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Abstract. An abstract should be given 



Quasi-periodic oscillations (QPOs) have been observed in 
accretion disks around neutron star, black hole, and white 
dwarf binaries with frequencies ranging from a few 0.1 Hz 
up to 1300 Hz. Recently, a correlation between their low- and 
high-frequency components was discovered and fitted with a 
single law, irrespective of the nature of the compact object. 
That such a relation holds over 6 orders of magnitude strongly 
supports the idea that the physical mechanism responsible for 
these oscillations should be the same in all binary systems. 

We propose a new model for these QPOs based on forced 
oscillations induced in the accretion disk due to the stellar 
magnetic field. First, it is shown that a magnetized accretion 
disk evolving in a rotating nonaxisymmetric magnetic field an- 
chored to a neutron star will be subject to three kinds of reso- 
nances: a corotation resonance, a Lindblad resonance due to a 
driving force, and a parametric resonance due to the time vary- 
ing epicyclic frequencies. The asymmetric part of the field is 
assumed to contain only one azimuthal mode m > 1. We fo- 
cus on the m = 1 disturbance, which is well studied for an 
inclined dipolar rotator; but our results are general and easily 
extend to m > 1 . However, the radial location of the resonances 
will be affected by this number m. For instance, with an m = 1 
asymmetric structure, the resonances reach regions very close 
to the innermost stable circular orbit (ISCO) and can account 
for observations of kHz-QPOs at frequencies as high as 1200- 
1300 Hz. If we replace the dipolar by a higher order multipolar 
component, the resonance location is shifted to larger radius, 
implying lower QPO frequencies. To compare the MHD situa- 
tion with the hydrodynamical case, we also consider an m = 2 
component in the magnetic perturbation in order to prove that, 
at least in the linear regime, the conclusions in both cases are 
the same. In the second part of the paper, we focus on the linear 
response of a thin accretion disk, developing the density per- 
turbation as the sum of free wave solutions and non-wavelike 
disturbances. In the last part, we show results of 2D numerical 
simulations of a simplified version of the accretion disk con- 
sisting of a column of plasma threaded by a vertical magnetic 
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field. These simulations are performed for the Newtonian grav- 
itational potential, as well as for a pseudo-general relativistic 
potential, which enables us to explore the behavior of the reso- 
nances around both rotating neutron stars and black holes. We 
found that the density perturbations are only significant in the 
region located close to the inner edge of the disk near the ISCO 
where the magnetic perturbation is maximal. They induce fluc- 
tuations in the density which persist over the whole time of 
the simulations and are closely related to the spin of the mag- 
netic perturbation. It is argued that the nearly periodic motion 
induced in the disk will produce high quality factor QPOs. 

Key words. Accretion, accretion disks - MHD - Instabilities 
- Methods: analytical - Methods: numerical - Stars: neutron 



1. INTRODUCTION 

To date, more than 20 accreting neutron stars in Low Mass X- 
Ray Binaries (LMXBs) are known to exhibit rapid variability 
in their X-ray fluxes. These high frequency quasi-periodic os- 
cillations (QPOs) show very strong similarities in their shape, 
as well as in their amplitude, and span frequencies from 300 Hz 
to about 1300 Hz (see van der Klis 2000 for a review). 

The comparison between high frequency QPOs (HFQPOs) 
in black holes and neutron stars shows a scaling proportional 
to 1/M*, where M* is the mass of the compact object, suggest- 
ing that these oscillations are due to the orbital motion in the 
accretion disk near the innermost stable circular orbit (ISCO) 
(McClintock & Remillard 2003). Moreover, the QPO frequen- 
cies depend on the spectral state of the source and are corre- 
lated with the accretion rate. QPOs are observed in the tail of 
the power law spectrum indicating that they are not due to ther- 
mal emission but probably related to Comptonization of soft 
photons. The saturation of the spectral index was explained by 
Titarchuk & Fiorito (2004) in the context of coronal heating 
and phase transition. 

Therefore, any QPO model must take not only oscilla- 
tions in the accretion disk into account but also the way 
photons propagate to the observer, including photon scatter- 
ing and general relativistic effects such as beaming, Doppler 
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shift, and aberration. This was investigated by Schnittman 
& Bertschinger (2004), who used exact integration of the 
geodesic motion of hot spots in the Kerr spacetime and looked 
at the harmonic content of the signal received by an observer 
in the asymptotically flat spacetime. 

In black hole binaries, the HFQPOs appear in pairs, as in 
the neutron star system, but their frequency remains uncorre- 
cted to the accretion rate. They act as a voiceprint at a ratio 3:2 
that favors a resonance mechanism. For example a parametric 
resonance between vertical and horizontal epicyclic frequen- 
cies was suggested by Kluzniak et al. (2004). Another mecha- 
nism proposed by Lee et al (2004) is resonance due to an ex- 
ternal driven force. 

The physics of accretion onto a magnetized neutron star has 
been described by Ghosh & Lamb (1978, 1979a, 1979b). In 
this picture close to the star's surface, the magnetic field chan- 
nels matter to the polar caps and exerts a torque on it, which 
accelerates or decelerates the neutron star depending on the an- 
gular momentum exchange rate. This is the standard model of 
magnetized accretion onto a neutron star. However, in order to 
explain pulsating X-ray sources, Anzer & Borner (1980) pro- 
pose an alternative view in which the magnetic moment of the 
neutron star and the accretion disk rotation axis are perpendic- 
ular. They also studied the influence of the Kelvin-Helmholtz 
instability on such a configuration (Anzer & Borner 1983). 

When the magnetic moment is misaligned with respect to 
the spin axis of the star, Vietri & Stella (1998) show that some 
inhomogeneities, treated as diamagnetic blobs, could be lifted 
off the disk and suffer the Lense-Thirring precession. 

Lai (1999) has identified a magnetically driven warping 
and precession instability that occurs in the inner region of 
an accretion disk. This occurs by interaction between the sur- 
face current density in the disk with the stellar magnetic field 
and produces motions at frequencies much lower than the or- 
bital one. It can, therefore, account for the low frequency 
QPOs (LFQPOs), as does the usual Lense-Thirring precession 
(Stella & Vietri 1998). Magnetic precession, combined with 
the relativistic dragging effect, was applied to weakly mag- 
netized neutron stars to explain these LFQPOS in LMXBs, 
(Shirakawa & Lai 2002a). The study of the global magnetic 
warping/precession process for strongly magnetized neutron 
stars can also explain the millihertz QPOs in some systems, see 
for example Shirakawa & Lai (2002b). The Rayleigh-Taylor in- 
stability associated with Rossby waves and the rotational split- 
ting gives another self-consistent description of the QPO phe- 
nomenon (Titarchuk 2003). 

It was realized by Kato (2001a, 2001b) that non- 
axisymmetric g-mode oscillations can be trapped in a thin rel- 
ativistic disk. These are excited by the corotation resonance. 
The effect of non-linear coupling between oscillations in the 
disk and warp is examined by Kato (2004). However, his study 
does not include the rotation of the compact object. 

Recent observations of accretion disks orbiting white 
dwarfs, neutron stars, or black holes have shown a strong cor- 
relation between the low and high frequency QPOs, LFQPOs 
and HFQPOs respectively, (Psaltis et al. 1999, Mauche 2002, 
Belloni et al 2002). This relation holds over more than 6 or- 
ders of magnitude and strongly supports the idea that the QPO 



phenomenon is a universal physical process independent of 
the nature of the compact object. Indeed, the presence or ab- 
sence of a solid surface, a magnetic field, or an event hori- 
zon plays no relevant role in the production of X-ray vari- 
ability (Wijnands 2001). This correlation has instead been ex- 
plained in terms of the centrifugal barrier model of Titarchuk 
et al. (2002) where the LFQPOs are associated with magneto- 
acoustic waves, and the HFQPOs correspond to the Keplerian 
motion in the disk. 

In this paper we propose a new resonance in the star-disk 
system arising from the response of the accretion disk to a non- 
axisymmetric rotating magnetic field. The paper is organized 
as follows. In Sec. 2, we describe the initial equilibrium state 
where the accretion disk is embedded in an axisymmetric mag- 
netic field anchored to the neutron star. Furthermore, we as- 
sume that the magnetic moment is aligned with the rotation 
axis of the star. We then estimate the field perturbation induced 
by a small inclination of the magnetic moment. In Sec. 3, we 
look for the resonance conditions and elucidate the meaning 
of the resonances. For a thin and weakly magnetized accretion 
disk, assuming no warp and only motion in the orbital plane 
of the disk, we investigate the precise linear behavior of the 
disk which is presented in Sec. 3.2. In Sec. 5, we perform 2D 
numerical MHD simulations of the disk and show that some 
resonances persist on a long timescale with a narrow radial 
extension. This is done for Newtonian, as well as for pseudo- 
general relativistic gravitational potential describing the Kerr 
geometry in a Newtonian way adapted to a rotating neutron 
star. Conclusions are drawn in Sec. 6. 

2. THE INITIAL CONDITIONS 

2. 1 . The stationary state 

In the equilibrium state, the accretion disk evolves in a perfectly 
spherically symmetric gravitational field and an axisymmetric 
dipolar magnetic field, both due to the compact object. We as- 
sume that the magnetic moment of the star and its rotational 
axis are aligned. The disk is assumed to lie in the equatorial 
plane of the star so that it experiences an axisymmetric station- 
ary field. By adopting a cylindrical coordinate system labeled 
by (r,<p,z), we can therefore describe all physical quantities 
with only a (r, z) dependence. Neglecting the matter inflow be- 
cause in the thin disk approximation v r <sc v v and because we 
only focus on oscillations around a stationary state, the den- 
sity p, pressure p, current density j, velocity v, and magnetic 
field B at equilibrium are thus given by : 



P = P(r,z) (1) 

P = P(r,z) (2) 

v = r£l(r,z)e v (3) 

B = B r (r,z)e r + B z (r,z)e z (4) 

j = jc P (r,z)e ip (5) 



The magnetic field has two different origins : 

- one component B* is attached to the neutron star and cor- 
responds to the rotating stellar dipolar magnetic field ; 
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- the second component B c i of the magnetic field is induced 
by the fluid motion in the accretion disk and follows from 
Ampere's law rot B d = p j. 

The total magnetic field is then simply the sum of these two 
components B = B t + Bd- 

In the radial direction, the gravitational attraction from the 
compact object g is balanced by the centrifugal force, the 
Lorentz force, and the pressure gradient, whereas in the ver- 
tical direction we simply have hydromagnetic equilibrium : 



K dp . n 3 ( B\ 

-p—=Pgr--Z-+]<pB z = pgr- — \p+ 

r or or \ 2p 



P8z 



Yz~ hBr 







(6) 
(7) 



To specify a given equilibrium state, we prescribe the initial 
density and current, or, equivalently, the magnetic field in the 
disk and the strength of the dipolar magnetic field of the neu- 
tron star. Assuming a thin accretion disk and a low /3-plasma 
parameter, defined as 



B 2 /2p 



(8) 



the pressure gradient and Lorentz forces will be negligible 
compared to the gravitational force. The motion will then only 
deviate slightly from Keplerian rotation. We then perturb the 
system by adding a small inclination angle \ between the mag- 
netic and rotational axis of the star. As a consequence, and due 
to the magnetic perturbation, the disk will leave its equilibrium 
state and start to oscillate around the equatorial plane. We next 
give the expression for the stellar magnetic field perturbation. 

2.2. Perturbed dipolar magnetic field of a rotating star 

In the case of an aligned rotator,^- = 0, the magnetic field is ax- 
isymmetric. Therefore, in the disk plane, there is only a vertical 
component B z . When the magnetic moment p is inclined by an 
angle x + relative to rotation axis Q„ an additional magnetic 
component dragged by the star's rotation arises. This perturba- 
tion is responsible for the oscillation in the disk. We derive its 
expression below. The magnetic moment can be developed into 
a parallel (along e z ) and a perpendicular (along e*) component 
such that 



p = p sin^e, + p cos^c z . 



(9) 



where e* is a radial unitary vector corotating with the star. In the 
Cartesian coordinates of a distant observer at rest (e x ,e y ,e z ), it 
can be expressed as 



e* = cos(Q* f) e x + sin(O t f) e y 



(10) 



assuming that at t = e* = e x . We adopt a particular geomet- 
ric configuration corresponding to the situation where the disk 
is contained in the star's equatorial plane. We could general- 
ize to the case where magnetic moment, disk, and star rotation 
axis are all misaligned. However, the rotational frequency of 
the magnetic perturbation, which is its most important prop- 
erty, would not be affected. 



The angular rotation is assumed to be equal to £1* = Q* e z . 
The stellar dipolar magnetic field has the usual expression de- 
rived from the vector potential 



_ p fi A r 



(11) 



4-n r 3 

Therefore, the perturbation is 5A* - j^p sin^-^^. It rep- 
resents a perpendicular rotator with magnetic moment p ± = 
p sinx- Expressed in the cylindrical coordinate system, the 
components of the perturbation of the magnetic field are given 
by: 

3r 2 



6b: 



Pq P smx 



4 TV (r 2 +z 2 ) 3/2 



r 2 + z 2 



1 



cos(<^ - £2, f) 



6Bt = 



5Bi = ^~ 



4n (r 2 +z 2 ) 3/2 
Po P sin^ 



sin(<p - O, t) 
3rz 



cos(^) - Q, t) 



(12) 



(13) 



(14) 



47r (r 2 +z 2 ) 3/2 r 2 +z 2 
We see that in the equatorial plane the perturbation is 5B z t = 
0. Each component has a time function that depends 
on {cos / sm}(<p - Q« t). The nature of the driving is therefore 
an m — 1 azimuthal mode with an excitation frequency at 
a rate Q.„. This has to be contrasted with the m — 2 mode 
and, consequently, excitation frequency of 2£2» produced by 
the quadrupolar gravitational field (the lowest order perturba- 
tion in gravity for a hydrodynamical accretion disk for exam- 
ple). This situation will be studied in detail in a forthcoming 
paper. The nature of the perturbation has some consequences 
for the derivation of the resonance conditions, as shown in the 
following Sec. 3. 

3. LINEAR ANALYSIS 

The distorted magnetic field described above will kick the disk 
out of its equilibrium state. We study its linear response by 
treating the asymmetric part of the magnetic field as a small 
perturbation of the equilibrium state prescribed in Sec. 2. We 
start with the ideal MHD equations of the accretion disk with 
adiabatic motion given by : 



^+div(pv) = 
ot 

— +(v-V)v = g-—L +J AB 
ot p 



Dt\ p y 



o 



dB 

— = rot(vAB) 
rotB = p j. 



(15) 
(16) 

(17) 

(18) 
(19) 



All quantities have their usual meaning: p density of mass in 
the disk, v the velocity of the disk, p the gaseous pressure, y 
the adiabatic index, g the spherically symmetric gravitational 
field, and B the total magnetic field. 

3.1. Lagrangian displacement 

We then perturb the equilibrium state from Sec. 2 by introduc- 
ing the Lagrangian displacement Expanding to first order we 
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get for the Eulerian perturbations of the magnetohydrodynam- 
ical quantities attached to the disk : 



Sp = -div(p£) 

Sv = -^+v-V£-£-Vv 

dt b b 

Sp = -ypdivg-g-Vp 

SB = B V£-Bdiv£-£ VB 



(20) 
(21) 
(22) 
(23) 



—1 |- (r(B„ + Q V + SB?)) - —(B r + Q r + SB r j\ 



dSB[ dSBi 



dz 



dr 
d5B 



r \dr dcp 
while for the azimuthal component we get : 



(25) 



For a detailed derivation of these perturbed quantities and their 
associated boundary conditions as well as for the equations de- 
scribed in this section, we refer the reader to Appendix A. By 
also making allowance for a perturbation in the neutron star 
magnetic field and following the Frieman-Rotenberg analy- 
sis (Frieman & Rotenberg 1960), the Lagrangian displacement 
satisfies a second order linear partial differential equation given 
by : 

p^f-div ^(vV)v]-Vn-lB-V6 
-Lq VB + div(p£)(g + Sg) - P Sg = 
— [rot (B + Q + SB t ) A SB, + rot SB t A (B + Qj] (24) 

HO 

We have introduced the convective derivative by D / Dt = 
d t + Q 8^. The perturbation in the magnetic field is represented 
by the vector Q = rot(^ A B). The scalar fl = ypdivg + 
£ • Vp - j^B ■ Q corresponds to the opposite of the total 
(gaseous+magnetic) pressure perturbation. Compared to the 
equation obtained by Frieman & Rotenberg (1960), we have 
an additional term in the Lorentz force induced by the stellar 
magnetic perturbation depicted by the expression enclosed by 
the brackets after the — factor, terms containing 6B t . We also 
allow for a degree of freedom in the choice of an arbitrary grav- 
itational perturbation induced by the disk 5g. However, for the 
rest of this paper, we only focus on the consequences of an 
asymmetry in the magnetic field, leaving the gravity perturba- 
tion for future study. 

We make a clear distinction between the stellar magnetic 
perturbation SB, which is known at each time and given by the 
rotating magnetic dipole Eq. (12)-(14) or any kind of multipo- 
lar fluctuation associated with the star and those caused by the 
current in the disk SB in response to the former perturbation. 

Expressing Eq. (24) in cylindrical coordinates for the radial 
component leads to : 

d 2 & _ dn 

p — 1 nil h 

F Dt 2 F Dt dr 

-^div(p£)+pr£-V(n 2 )- 
p dr 

dQ r By dQ r B^Qy 

D r — - — + — — — Z h 

dr r dip r 



1 

Ho 
dQ r 



1 



SBl 



t , dB r Q v dB r dB r 
B z — + Q r — + — — + Q z — 

dz dr r dip dz 

{ j z {B r + Q r + SBO - |- (B z + Q z + 6BI)) 



dQ, 



1 

P0 

dB 



° ^ +o n D ^ r 1811 
p -5F +2p ^-bl--r&p- 

_ dQ v B<p dQ v By Q r 

B r — 1 1 h 

dr r dip r 



Qip dBy Qy B r 



dz 



dr 

1 

y"0 



dip 



dz 



SB" { d i w s 

-f \y r {r{B, + Q, + 6Bt))- 

^(B r + Q r + dB r j\ - 

{1 d d 1 

-■^(B z + Q z + SBO - j z (B v + Q V + 6Bt)\ 



Br+Qr I d dSB[ 

r \dr dtp 

1 dSB; dSBt 



(B z + Q z )\- 

\r dip dz 

and finally for the vertical oscillations we have : 



PTJ - "37 " — 



Dt 2 



dn 

dz 



1 

14) 



dQ B^ dQ dQ 

D r — — + — + D 7 — h 

dr r dip dz 



dB z Q^ dB 



dB, 



Qr — ~ — — - (A — 



dr 



r dip 
1 

Ho 



dz 



-div (pg))g z + 



SBt\-^-(B z + Q z + SBD- 
yr dip 

Q 

-(B^ + Q^ + SBt) 
- (B r +Q r + SBl) - - (B z + Q z + SB? 



1 1 dSBl 
(8, + C ,)(f^ 



dSB 
dSBV 



+ 



(26) 



(27) 



We emphasize the fact that these expressions contain terms of 
the form ^ <5B», = r,ip,z) that are second order with re- 
spect to the perturbation and that, therefore, could be neglected 
as long as SB* remains weak. But in doing so, we suppress the 
parametric resonance to be studied in more detail below. This 
resonance becomes relevant when the perturbed stellar mag- 
netic field is on the same order of magnitude as the aligned 
component. Depending on the magnitude of the perturbation, 
which we suppose small compared to the background field, 
this resonance will develop on a timescale closely related to 
the amplitude of the perturbation so should not be ignored. The 
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terms SB' t 8B{ are also second order when the magnetic asym- 
metry is small. 

Finding an analytical stability criterion for this system is a 
complicated or perhaps even an impossible task. Furthermore, 
we cannot apply the classical expansion in plane wave solutions 
leading to an eigenvalue problem. The presence of some coef- 
ficients varying periodically in time, like SB*, prevents such a 
treatment. However, the problem can be cast in a more conve- 
nient form if we decouple the oscillations in the orbital plane 
from those perpendicular to it, as is done in the following sub- 
section. 

3.2. Detailed linear analysis of the thin disk 

If we assume that there is no warping in the disk and ne- 
glect the gradient and displacement in the vertical direction, 
we can perform a detailed and accurate analysis of the mo- 
tions in the disk. Indeed, by setting £ z = and d/dz = in 
Eq. (25) and (26), the linearized MHD equations become much 
more tractable. Moreover, we first assume a simple geometric 
structure in which the disk is replaced by a column of plasma 
threaded by a vertical magnetic field and then allow only mag- 
netic perturbation in this direction. This situation seems far 
from the real system of an oblique rotator; nevertheless, it gives 
a good first insight into the oscillations and resonances arising 
there. 

As a result, setting B r = B v = SB r t = SBt = we find for 
the magnetic perturbation 



Qr 

Qz 



d_ 

dcp ' 



(rB z f r ) + ^(Bz&) 



and for the opposite of the total pressure perturbation : 

dp B z Q z 



n 



y_p_ 

r 



dr (r&) + ^ 



+ 6 



dr 



pa 



(28) 
(29) 

(30) 



(31) 



We seek solutions by writing each perturbation, such as the 
components of the perturbed magnetic field SB (i.e. those of 
Q), those of the Lagrangian displacement those of the per- 
turbed velocity Sv, the perturbed density Sp, and the perturbed 
pressure Sp as 



X(r, <p, t) = Re[X(r) e Kn " f -°'' ) ~i, 



(32) 



where m is the azimuthal wavenumber and cr the eigenfre- 
quency of the perturbation related to the speed pattern Q. p by 



cr = mQ. p . Therefore 
1 d 



~ i o i ~\im 
Qz = — T (rB z { r )-—B z ^ 
r dr v ' r 



n = P c 



2 

tiuiz 



1 d 



1 U I ~ \ I III ~ 

-ro-M) + T^ 



, d 



B 1 



+ & 3-r[? + 2p- 



(33) 



(34) 



c 2 +cl z is the vertical fast magneto-acoustic wave speed. 



vertical magnetic field are respectively : 



cl = 



y_p_ 
p 

Pop 



(35) 
(36) 



Neglecting quadratic expressions in the excitation term, 
like Q Z SB\, the azimuthal Lagrangian displacement Eq. (26) 
is solved by : 



%tp — 2 



m B Z 5B\ 



r P Po 



+ l2 ^ + -pdr[ P + ^-oP &+ 



mc maz 9 



(37) 



Due to the differential rotation of the disk, the eigenfre- 
quency cr appears Doppler-shifted. It is therefore convenient 
to introduce the frequency w and a second one w» taking the 
speed of the magneto-sonic wave into account by 



a> = cr - mil 

2 2 

2 2 m C maz 

or = or — . 



(38) 
(39) 



Strictly speaking, corotation is achieved when co„ = 0, such that 
there exist two solutions. When far from the corotation points, 
because of the weakly magnetized thin accretion disk approx- 
imation 0(c 2 s ) = 0(cl z ) <K 0(r 2 cj 2 ) and therefore co 2 . « oj 2 . 
Substituting the expressions for Q z , Eq. (33), for fl, Eq. (34) 
and for f^, Eq. (37), in the radial equation Eq. (25), we obtain 
to lowest order : 



d 2 l r dln(rpc 2 maz ) d\r 



dr 2 



dr 



dr 



A B z dSBl 

' POP 



dr 
(40) 



The radial epicyclic frequency of a single particle as measured 
by an observer in the rest frame is 



K 2 =A£l 2 + r — 
dr 



(41) 



(see Appendix A). The solutions of this equation consist of free 
waves corresponding to density perturbations propagating in 
the disk without magnetic asymmetry related to the homoge- 
neous part and of non-wavelike disturbances due to this asym- 
metric component related to the inhomogeneous part. 

3.2.1 . Free wave solutions 

When looking for free wave solutions to Eq.(40), we can 
set SBl = 0. Thus 



d 2 l | dln(rpcj az ) dl 
dr 2 dr dr 



+ {a? -*?)| r = 0. 



(42) 



This is the generalization to MHD of the sound wave that prop- 
agates in the hydrodynamical disk. We just need to replace 



The sound speed and the Alfven velocity associated with the the sound speed c 2 by the fast magneto-acoustic speed c„ 
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Introducing a new unknown if/ = | r ^Jr p c 2 maz , it satisfies a kind 
of Schrodinger equation 

i//"(r) + V(r)if/(r) = (43) 



where the potential is given by : 
V(r) 



2 2 
OJ - KZ 



(44) 



A first guess for free wave solutions is given by the WKB ex- 
pansion as follows : 



(45) 



Putting this approximation into Eq. (43), the dispersion relation 
for the fast magneto-acoustic waves propagating in the orbital 
plane of the accretion disk is given by : 

. .2 _ ..2 , 2 i 2 



k 2 + c 2 k L 

K r + c maz K 



(46) 



Free waves can only propagate in regions where <x> - k 2 = 
V( r ) c maz — 0- The boundary between the propagating and 
damping zones is defined by the inner and outer Lindblad 



radius r 



in /our 



defined by V(r™ lout ) 



Appendix B, we can find a better approximation to the solu- 
tion of Eq. (43) which is valid even for r « r ™l° ut _ p or th e inner 
Lindblad resonance of interest here, we introduce the following 
function oj\, writing r L -r": 



a>i(r) 




for r < ri 



for r > r^. 



(47) 



(48) 



The function iff is then a linear combination of the 2 linearly 
independent solutions: 



Mr) 



Mr) = 



Ai{oJi(r)) 

Bijmjr)) 
yfW)\' 



(49) 



(50) 



Furthermore, we require the solution to remain bounded, which 
leads to C2 = 0. Thus the solution for the Lagrangian dis- 
placement is % r = C\ iffi(r)/ ^rpc 2 mav At the inner boundary 
of the accretion disk, in accordance with the simulations per- 
formed later on, the density perturbation should vanish. This 
is expressed as dp = 0. To the lowest order consistent with 
our approximation, the Lagrangian radial displacement £ r must 
satisfy 



d\n(rp) 



dr 



n 

+ 2m 



URi) = 0. 



(51) 



This last condition determines the eigenfrequencies u as a 
function of the azimuthal mode m. For any m, there is an infi- 
nite set of eigenvalues. However, the corresponding eigenfunc- 
tions become more and more oscillatory, implying larger and 
larger wave numbers. In the numerical applications, we restrict 
our attention to the first few that also correspond to the highest 
possible eigenvalues cr. 



3.2.2. Newtonian vs general-relativistic disk 

We emphasize that the qualitative results presented in this pa- 
per are independent of the nature of the spacetime, be it flat as 
in the Newtonian gravitational potential or be it curved as in 
the Kerr metric. We make this point clear by showing results 
according to both spacetime structures. 

Let's start by recalling some facts about accretion disks or- 
biting in the Kerr metric. When the inner edge of the accretion 
disk reaches a few gravitational radii, general relativistic ef- 
fects become important. The degeneracy between the three fre- 
quencies is lifted, namely the orbital £2, the radial epicyclic K r , 
and the vertical epicyclic k z frequencies. Their value depends 
not only on the stellar mass M, but also on the angular mo- 
mentum of the star a„. Indeed we distinguish 3 characteristic 
frequencies in the accretion disk around a Kerr black hole (or 
equivalently a rotating neutron star): 



the orbital angular velocity: 
0(r, a,) 



0. Using the results of _ 



r 3 ' 2 + a t 
the radial epicyclic frequency: 



16 a« a 2 
K r {r, a.) = &{r, a») -W 1 - - + 8 - 3 



- the vertical epicyclic frequency: 



K z (r, a,) = Q.(r, a,) -W 1 - 4 — + 3 ^ 



(52) 



(53) 



(54) 



The parameter a, corresponds to the angular momentum of the 
star in geometrized units. For a neutron star, it is given by a„ = 



GMl 



a. 

We have the following ordering : 



Q > k z > K r for a, > 
k z > Q,> K r for a, < 

3.2.3. Results 



(55) 
(56) 



The eigenvalues for the density waves are shown with decreas- 
ing magnitude in Table (1). This holds for a neutron star with 
angular velocity v, = Q.*/2n = 100 Hz. We compared the 
Newtonian case with the Schwarzschild metric. The highest 
speed pattern, given by cr/m, never exceeds the orbital fre- 
quency at the ISCO. Note that in the m = 1 case, the inner 
Lindblad resonance does not exist for the Newtonian poten- 
tial, and only a few free wave solutions can be found for the 
Schwarzschild geometry. 

Some examples of the corresponding eigenfunctions for 
density waves are shown in Fig. 1 with arbitrary normalization. 
Each of them possesses its own inner Lindblad radius depend- 
ing on the eigenvalue. They are distinguished by the number of 
radial nodes they possess, starting from no node at all, which 
corresponds to the highest speed pattern. 

Because the pattern speed of the first spiral waves is close 
to the ISCO, we suggest that these free wave solutions can be 
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Table 1. The first eight largest eigenvalues a for the free wave so- 
lutions of Eq. (43). Values are normalized to the frequency of the 
ISCO, £lj sc0 = 6~ 3/2 . Results are given for 3 azimuthal modes m = 
1,2,5 for the Newtonian, as well as for the Schwarzschild gravita- 
tional field. Where no inner Lindblad resonance exist, this is repre- 
sented by - . 



Eigenvalues cr/fl,, 



Newtonian 


Schwarzschild 


m = 1 m = 2 


m = 5 


m = 1 


m = 2 


m = 5 


0.745 


3.286 


0.3727 


1.139 


3.643 


0.484 


2.564 


0.126 


0.704 


2.802 


0.358 


2.155 


0.0528 


0.504 


2.330 


0.278 


1.870 


0.0217 


0.382 


2.008 


0.220 


1.652 


0.010 


0.298 


1.765 


0.177 


1.476 




0.236 


1.571 


0.143 


1.330 




0.189 


1.411 


0.116 


1.206 




0.152 


1.276 





Fig. 1. Density wave perturbation in the disk caused by the free wave 
propagation for the azimuthal mode m = 2. Some examples are shown 
for different eigenvalues and for Newtonian geometry, on the left, as 
well as for Schwarzschild geometry, on the right. The vertical bar indi- 
cates the location of the inner Lindblad resonance. The normalization 
of the eigenfunctions is arbitrary. 



due to the varying accretion rate. For instance when the accre- 
tion process is enhanced, the inner edge moves closer to the 
ISCO. As a result, the highest eigenvalue of the free waves also 
increases and therefore v 2 shifts to higher frequencies as well 
as Vi. 

Examples of this mechanism are shown in Fig. 2. We es- 
timate the highest eigenvalue for the modes m = 2, 5 in the 
Newtonian, as well as in the Schwarzschild gravitational po- 
tential. The pattern speed of the wave solution D. p , normalized 
to the frequency of the ISCO, Q.j SC0 , is shown as a function of 
the location of the inner boundary of the accretion disk /?,„. 
In each case, we observe an increase in the eigenvalue when 
the disk's inner edge moves closer to the compact object, cor- 
responding to an increased accretion rate. When the ISCO is 
reached, the boundary condition does not change anymore and 
the eigenvalue saturates to its final value. Thus the kHz-QPO 
frequencies associated with the speed pattern of the wave satu- 
rate. This has been observed in some LMXBs, as reported for 
instance in a paper by Zhang et al. (1998). 




Fig. 2. Variation of the highest eigenvalue, corresponding to the eigen- 
function having no node, as a function of the location of the inner edge 
of the disk R jn . There is a monotonic increase as the disk approaches 
the ISCO at R jn = 6R g . Results are shown for the m = 2, 5 modes in the 
Newtonian (N) and Schwarzschild (S) spacetime, red and blue curves 
respectively. The gravitational radius is defined by R g = GM,/c 2 . 



associated with the highest kHz QPO More precisely, there ex- 
ists one eigenfunction which possesses no root, namely the one 
with the highest eigenvalue. Eigenfunctions with several nodes 
in the radial distribution will emit no significant radiation since 
on average the fluctuations will cancel out. Therefore, we ex- 
pect to see only density fluctuations with no or very few radial 
nodes. The most interesting candidate is the one with no node, 
a situation that also corresponds to the highest eigenvalue. 

Moreover, while propagating in the perturbed rotating grav- 
itational field of the star, this spiral wave experiences a sinu- 
soidally varying gravity at a rate O p -Q», where Cl p is the speed 
pattern of the wave. This induces a modulation in the shape of 
the wave which will be reflected in the power spectrum density 
as a kind of beat phenomenon. As a result, in our model, we 
associate the highest kHz-QPO to v 2 = Q. p and the lowest kHz 
QPO to the beat frequency v\ = Q p - £2,. The peak separation 
is then Av = £2„. 

In a real accretion disk, the precise location of the inner 
edge does not necessarily reach the ISCO. The magnetic field 
or the radiation pressure can prevent the disk from doing so. 
Indeed, the boundary of the inner part of the disk can fluctuate 



3.2.4. Non-wavelike perturbation 

We now solve the full inhomogeneous Eq. (40) in order to find 
the solution corresponding to the non-wavelike perturbation. 
In this case, the eigenvalue is known: cr = mO«. We have to 
solve the second order ordinary differential equation for if/ with 
the appropriate boundary conditions Eq. (51). Following the 
results given in Appendix B, the most general solution is: 

if/ r (r) = Ci iAi(r) + C 2 1A2O) + 7r signer)) x 

f ^ l (r)ifr 2 (s)- l / /l (s)^ 2 (r))F(s)ds. (57) 
•Jr L 

The constant Co is chosen such that the solution remains 
bounded for r » re 

C2 — lim n sign(<j/, ) I iffi(s) F(s) ds. (58) 

This integral is convergent because the function ip\ is expo- 
nentially decreasing with the radius. We can directly integrate 
this equation to obtain the shape of the density perturbation 
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between the inner edge of the accretion disk and the inner 
Lindblad resonance. The results for a particular choice of the 
density and magnetic profiles (discussed in Sec. 5) are shown in 
Fig. 3. This can be used later to check the numerical algorithm. 
It also gives a good estimate of the amplitude of these fluctua- 
tions because it is directly related to the strength of the stellar 
magnetic field asymmetry. The maximum density fluctuation is 
about 0.1%. The oscillations become evanescent upon reaching 
the inner Lindblad resonance. 



u ,.,iu-i 




W 11 f t% l-.l jn tl J 



Fig. 3. Non- wavelike disturbance in the inner part of the accretion disk 
due to the asymmetric component of the magnetic field. The potential 
is Newtonian and m = 2. The numerical values used are described in 
Sec.5. 



3.2.5. Counter-rotating disk 

The counterrotating disk can be solved in exactly the same way 
by integrating Eq. (40) to find the solution corresponding to the 
non-wavelike perturbation. The amplitude of the oscillations 
is shown in Fig. 4. The Lindblad resonances are no longer in 
the computational domain. An oscillatory structure appears in 
the whole disk. The maximum amplitude is about 1% of the 
average density. This can later be compared with the numerical 
simulations. 




Fig. 4. Non- wavelike disturbance in the inner part of the accretion disk 
due to the asymmetric component of the magnetic field. The potential 
is Newtonian, m = 2 and the disk is counterrotating. 



4. SIMPLIFIED ANALYSIS 

To get more insight into the nature of the resonances, we now 
focus on the displacement of the disk either in the equatorial 
plane or perpendicular to it. This means that we neglect the 
coupling between the oscillations occurring in perpendicular 
directions. This simplification will help us to understand the 
origin of the instability and to derive some resonance condi- 
tions. 

To study the motion in the vertical direction, set- 
ting (£ r , i; v ) = we find: 

D 2 & d 17 2 B r SB r A 

+ 



p 



Dt 2 



d 

dz 



d_ 

d~z 



5B r , 2 + SBV B r 5B r „ 



y"0 



dz \ fi Q dz * z ) 



(59) 



See Appendix A for the details. The Alfven velocity associ- 
ated with the radial magnetic field and the radial fast magneto- 
acoustic wave speed are, respectively: 

i B 2 



c nr = 



Pop 

2 



(60) 
(61) 



The same can be done for the orbital motion by setting £■ = 
0, which leads to a similar expression : 



Dt 2 



d_ 
~dr 



P C maz 

d (SB? 
dr\2^o Po 



P0 

B Z 6B 



r or 



l\ d (6BI dB, \ 



Eq. (59) and (62) look very similar, the discrepancy coming 
only from the difference between the planar and the cylindrical 
geometry (terms containing r). In the two first terms of Eq. (59) 
and (62) we recognize an MHD wave propagation in a tube of 
spatially varying but time independent cross section (Morse & 
Feshbach 1953). The first and third terms put together form an 
harmonic oscillator at the epicyclic frequency. So the three first 
terms are a generalization of the Klein-Gordon equation and do 
not give rise to any kind of instability. The interesting parts are 
those containing the perturbation in the magnetic field SB,. 

Further, we neglect the magneto-acoustic wave propaga- 
tion because of the weakly magnetized thin disk approxima- 
tion, and also assume a homogeneous medium for which the 
vertical Lagrangian displacement and density are independent 
of the altitude z in the first case and where the density and ra- 
dial Lagrangian displacement do not depend on the radius r in 
the second case. Note also that terms like 6B' t are second order 
with respect to the perturbation and should be neglected. These 
simplifications make the equations much more easier to manip- 
ulate without losing the most interesting physical aspect con- 
tained in the resonance terms. We can then simplified Eq. (59) 
and Eq. (62) to obtain : 



D 2 £ 
Dt 2 



2 _d_l^K_ d^- 

z dz\HoP dz 
2 d 16B% 8B Z 
r dr \// p dr 



d lB r 5B[ 
dz\ flop 
d IB Z 8B\ 
dr\ flop 



(63) 
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(64) 

Recall that the magnetic perturbation is dragged with the star's 
rotation. Therefore, when looking at a given fixed point in 
the accretion disk, the coefficients of the equations contain- 
ing 6B' t will vary periodically. Thus, we recognize a kind of 
Hill equation corresponding to an oscillator with periodically 
time-varying eigenfrequency and subject to a periodic driven 
force. It is well known that this type of equation shows what 
is called a parametric resonance. Moreover, due to the rotation 
of the star, this perturbation will vary sinusoidally in time. In 
fact, in the frame locally corotating with the disk, (convective 
derivative), the modulation follows a harmonic temporal be- 
havior such as : 



cos(ip - (Q. - Q„) f) , sin(<£> - (Q - O t ) f). 



(65) 



The Hill equation specializes to Mathieu's equation for which 
we can derive the resonance conditions analytically. The solu- 
tions to Mathieu's equation written in the form y"(t) + w^(l + 
h cosyf)y(f) = will indeed become unbounded if y = 2^ 
where n is an integer, (Landau & Lifshitz 1982). 

The equation satisfied by each component of the 
Lagrangian displacement when looking at a given azimuth <po 
is therefore in the form : 



r (o + 



/ 



{co"}^ ' 



f(0 = 
(O - Q„) 0, 



(66) 



where f{5B'^) and h{5B'^) are two constants depending on the 
strength of the stellar magnetic perturbation. 

From this analysis we expect three kinds of resonances cor- 
responding to, (see Appendix C for a detailed analysis) : 

- a corotation resonance at the radius where the angular ve- 
locity of the disk equals the rotation speed of the star. This 
is only possible for prograde disks. The resonance condi- 
tion to determine the corotating radius is simply 



Q = Q, 



(67) 



- an inner and an outer Lindblad resonance at the radii where 
the radial or vertical epicyclic frequency equals the fre- 
quency of the dipolar magnetic perturbation as measured 
in the locally comoving frame. Due to the m = 1 structure 
of this perturbation we find the resonance condition to be 



7z- 



(68) 



- a parametric resonance related to the time-varying 
epicyclic frequency, (Hill equation). The rotation of the 
star induces a sinusoidally time varying epicyclic frequency 
leading to the well-known Mathieu equation. The reso- 
nance condition is then 



in- n.i = 2^. 

n 



(69) 



Note that the driven resonance occurs at the same radius as 
does the parametric resonance for n = 2. However, their growth 
rates are different, because forcing causes a linear growth in 



time of the amplitude while parametric resonance causes an 
exponential growth. 

More generally, for a perturbation of arbitrary azimuthal 
mode m, the resonance conditions are given by : 

- for the inner and outer Lindblad resonances : m \Q. - Cl,\ = 

K r/Z 

- for the parametric resonance : m \Q - Cl,\ -2-^-. 

The conclusion of this study is that the locations of the reso- 
nances are determined solely by the rotation rate of the neu- 
tron star. Fluctuations of the accretion rate are not included in 
this picture. However, we will show by a more detailed linear 
analysis that there are some free-wave solutions for the density 
perturbation for which the eigenvalue can be correlated to the 
accretion rate in accordance with observations, see Sec. 3.2. 

4.1. Results 

4.1 .1 . Newtonian disk 

From the resonance conditions derived previously, Eq. (67)- 
(69), we can estimate the radii where each of these resonances 
occurs. Beginning with the Newtonian potential, it is well 
known that the angular velocity, the radial and epicyclic fre- 
quencies are all equal. For a weakly magnetized thin accre- 
tion disk, this statement remains true so that n « K r « k z . 
Distinguishing between the two signs of the absolute value, we 
get the following rotation rate for the parametric condition (69): 



n 



n 

n* n + 2 
Q n 

nT ~ n-2 



= 1/3,1/2,3/5,2/3,5/7,... 
= -1, -,3,2,5/3,... 



(70) 
(71) 



As a consequence, each resonance will occur in the frequency 
range [Q./3, 3 n.]. Assuming that the QPOs are related to the 
orbital motion and taking into account that the spatial structures 
of the resonances are formed by m = 1 patterns, we expect 
to see QPOs in the same range [v,/3, 3 v.] Hz where we have 
introduced the spin frequency by v* = n»/27r. For instance, 
for a 300 Hz spinning neutron star, we can expect QPOs be- 
tween [100,900] Hz. However, observations have shown that 
QPOs with frequency up to 1300 Hz exist for this kind of sys- 
tem. We will see in the next subsection that this difficulty can 
be solved by introducing the general relativistic effects of the 
innermost stable circular orbit. Table 2 gives an example of the 
expected orbital frequencies where the resonances are expected 
to occur for a neutron star spinning at 300 and 600 Hz. 

4.1.2. General relativistic disk 

The parametric resonance condition Eq. (69), split into the two 
cases depending on the sign of the absolute value, becomes: 



fl(r,a.)-2 f/ - v = n» 
n 

n(r, a,) + 2 = O* 

n 

For a typical neutron star, we choose: 



(72) 
(73) 
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Table 2. Value of the orbital frequencies at the parametric resonance 
radii for the first five orders n in the case of a Newtonian gravita- 
tional potential. The results are given for a 1 .4 M Q neutron star rotating 
at 300 and 600 Hz respectively. The value to the left of the symbol / 
corresponds to the absolute value sign taken to be - and on the right 

to be +. The sign indicates that no resonance condition can be 

satisfied. 



rank n Orbital frequency v(r, a,) (Hz) 





v, = 600 Hz 


v, = 300 Hz 


1 


-600 / 200 


-300 / 100 


2 


— / 300 


— / 150 


3 


1800/360 


900/ 180 


4 


1200/400 


600 / 200 


5 


1000/429 


500/214 



Table 3. Value of the orbital frequencies at the parametric resonance 
radii for the first five orders n in the general relativistic case. The re- 
sults are given for a 1.4 M Q neutron star rotating at 300 and 600 Hz 
respectively. The value to the left of the symbol / corresponds to the 
absolute value sign taken to be - and on the right to be +. 



rank n 


Orbital frequency v(r, a,) (Hz) 




Vertical 


Radial 




v, = 600 Hz 


v, = 300 Hz 


v, = 600 Hz 


v« = 300 Hz 


1 


— /200 


— / 100 


1596 / 220 


1332/ 106 


2 


— / 301 


— / 150 


1196/330 


800/ 159 


3 


1695 / 361 


885/ 180 


980/392 


573 / 190 


4 


1175/401 


597 / 200 


870/432 


480/210 


5 


988 / 430 


498/214 


808 / 459 


433 / 223 



- mass M, = 1.4M ; 

- angular velocity v» 

- moment of inertia /, 

- angular momentum a, 



0«/2tt = 300-600 Hz; 
= 10 38 kgm 2 ; 



G Ml 



The angular momentum is then given by a* = 5.79 * 10~ 5 Q*. 
For a given angular momentum a*, we have to solve equa- 
tions (72), (73). Solving them for the radius r and then deducing 
the orbital frequency at this radius, we get the results shown in 
Table 3. For the spin rate of the star we find a, = 0.109 - 0.218 
and so the vertical epicyclic frequency is close to the orbital 
one k z « £1 Thus for the vertical resonance, we are still close 
to the Newtonian case explored in the previous section. 

We have now overcome the problem faced in the 
Newtonian approximation. Indeed, even in the case of a 300 Hz 
spinning neutron star, we can expect frequencies as high as 
1332 Hz for the n = 1 rank, in accordance with observations. 
The general-relativistic effect caused by the presence of the 
ISCO is therefore a crucial characteristic that makes a clear 
distinction with the Newtonian potential. The ISCO has to be 
implemented in any realistic QPO model in order to get a good 
quantitative agreement with observations. Therefore, a fully 
general-relativistic simulation code should be developed in or- 
der to take this ISCO into account in a self-consistent manner. 
However, we point out that general relativistic effects are not 
directly responsible for the observed QPOs as claimed in other 
models (Lee et al 2004, Stella et al. 1999). They merely lead 



to quantitative improvements compared to the flat spacetime 
results. 



5. NUMERICAL SIMULATIONS 

The linear analysis described above is not the end of the story. 
It was shown that the disk becomes unstable at some preferred 
radii where the resonance conditions are satisfied. But this anal- 
ysis does not tell us what will happen to these resonances on a 
long time scale. Moreover, as in the hydrodynamical case, the 
non-linearities in the time evolution of the disk will give rise 
to a rich variety of physical processes. However, in this simu- 
lation, in order to check our numerical algorithm, we start with 
runs for which the perturbation remains weak, in order to keep 
a regime close to the linear one described in the previous sec- 
tion. The full non-linear effects arising from a strong magnetic 
field perturbation will be treated in a forthcoming paper. 

To achieve this goal, we performed 2D numerical simula- 
tions by solving the magnetohydrodynamical equations for the 
accretion disk Eq. (15)-(19). This is done by a pseudo-spectral 
: method for a simplified version of the 2D MHD accretion disk. 
-Indeed, an exact treatment of the magnetic field induced by the 
motions in the disk would involve solution of an elliptic partial 
differential equation related to Ampere's law. We don't want to 
go into such refinement in this paper but only emphasize the 
role of a rotating asymmetric magnetic field. Nevertheless the 
physics at hand remains the same. 

5.1. Boundary and initial conditions 

We expand every quantity in Fourier series in the azimuthal 
direction and in Tchebyshev series in the radial direction, 
whereas time integration is performed by an explicit fourth- 
order Runge-Kutta scheme. In this work we are not interested 
in the accretion process itself but focus only on the oscillations 
occurring in this system around its stationary state. We can 
then impose fixed boundary conditions in the radial direction 
by setting the radial velocity equal to zero at the two bound- 
aries. Furthermore, to avoid spurious reflections at the edges 
of the accretion disk, we force the radial velocity to decrease to 
zero in a narrow layer around these boundaries. We also require 
the density and the azimuthal velocity to reach their equilib- 
rium. This acts as an ad hoc absorption process damping out the 
strong oscillations which could arise in these regions. The non- 
reflecting boundaries conditions are therefore expressed as: 



p(r, tp) = f(r) (p(r, <p) - p e %r, tp)) + p e «(r, tp) 

v r (r,<p) = f(r)v r (r,<p) 

v<p(r,<p) = f(r) (v v (r, tp) - v^(r, tp)) + v e *(r, tp) 

B z (r, tp) = f(r) (B z (r, tp) - B e z \r, tp)) + B e z \r, tp) 



(74) 
(75) 
(76) 
(77) 



The quantities denoted by an X subscript correspond to the 
physical values before applying the non-reflecting boundary 
conditions. Those values at equilibrium are denoted by X eq . 
Function / acts as a filter, as its value is one everywhere ex- 
cept close to the inner and outer boundaries where it decreases 
to zero in a narrow transition layer. 
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We adopted the geometry of a cylinder of gas threaded by 
a vertical magnetic field and a radially directed gravitational 
field. The problem is invariant under translation in the direction 
of the magnetic field lines, which is chosen to be the z direction, 
d/dz = 0. Motions of the gas are only allowed in the equatorial 
plane (r, <p), v z = 0; therefore, we suppress the possibility of 
a warped disk and a possible precession of its orbital plane, 
leaving this to future work. The ideal MHD equations (15) and 
(16) then simplify to : 



-£ + - -jrirpvr) + - —(pvy) = 
at r Or r o<p 



dv r 



dv r v<p dv r 



1 dp j v B z 



dt + Vr dr + r dtp r ^ r p dr + / ) 



dv w 



—— + v r — — + — + 

at or r dtp r 



8f 



1 dp j r B z 
pr d(p p 



(78) 



(79) 



(80) 



Moreover, we adopt a polytropic equation of state for the 
plasma which is: 



p = K R y. 



(81) 



The induction equation (18) in this particular geometry trans- 
forms into an advection equation for the vertical component B z , 
similar to the mass conservation. This means that the magnetic 
flux is conserved during time evolution: 



dB 1 d 1 d 

— z - + - —(rB z v r ) + - —(B z v v ) = 
at r or r dtp 



dB r 

~dT 

dB v 

dt 



= 



= 0. 



(82) 
(83) 
(84) 



Inspecting these expressions, we conclude that there is no ra- 
dial or azimuthal magnetic field generated during the evolution 
of the disk. Because at the initial conditions B r = B^ = 0, 
these components will remain equal to zero during the whole 
simulation. Consequently the divergence free condition of B is 
automatically satisfied by: 



divB 



dBz 

dz 



0. 



(85) 



Finally, applying Ampere's law, Eq. (19), we get the current 
density in the disk by: 



Jr = 



J_dB, 
fx r dip 

__L ^1 

Jlp A'o dr 

h = o 



(86) 

(87) 
(88) 



We expand every quantity in a Fourier series in the az- 
imuthal direction and in a Tchebyshev series in the radial di- 
rection. Moreover, to avoid the spectral aliasing that arises 
from the non-linearities of the equations, we use the 2/3-rule 
to completely wash out the highest frequencies (Gottlieb and 
Orszag 1977, Canuto et al. 1990) (desaliasing process). 



We use geometrized units in which G = c = 1. 
Distance is measured in units of the gravitational radius given 
by R g = G M t /c 2 . Moreover, the simulations are done for a star 
with M, = 1, so that in the new units we have R g = 1. The star 
is assumed to be perfectly spherical with the three main axes 
all equal R x = R y — R z , and the perturbation comes from the 
magnetic field alone. The adiabatic constant is equal to y = 5/3 
and typical resolution is 256 x 32. 

We checked that this resolution is sufficient by running 
a simulation with twice as many points in both directions, 
namely 5 12 x 64. We found no changes in the solutions. Indeed 
inspecting the Fourier-Tchebyshev coefficients, only half of the 
radial coefficients and a few of the azimuthal coefficients in the 
2D expansions are significantly different from zero. Choosing 
a higher resolution will not lead to any improvement in the so- 
lution but only to excessive computation time. 

The initial density profile is prescribed as 



PoW 



10" 



(89) 



Initially the contribution to the magnetic field is given by the 
star and the disk as follows 



B* = -5 



io- 



B a 



io- 



while the stellar magnetic perturbation is given by 

5B\ = 5b B* cos(m (tp - Q* f)) 
8b = 0.001. 



(90) 
(91) 



(92) 
(93) 



The r-dependence of the magnetic perturbation Eq. (92) is 
chosen to mimic the rotation of a misaligned rotator rotat- 
ing at a rate Q.„ with a typical decay as r -3 . All other com- 
ponents are equal to zero. The azimuthal mode of the mag- 
netic field perturbation 5B* can be chosen arbitrarily. We use 
the case m — 2 to directly compare with the hydrodynamical 
simulations. However, the m = 1 case is better adapted to the 
true dipolar structure of the stellar magnetic field. We will give 
some examples for this situation, too. 

With these prescribed initial conditions, the thin disk ap- 
proximation is fairly good. Indeed the disk ratio H/R, H typical 
height and R radius of the disk, plotted for the Newtonian and 
Schwarzschild geometry in Fig. (5), never exceed 0.1 in any 
run. Moreover, the disk is threaded by a weak magnetic field 
except very close to the inner edge. The plasma B parameter is 
shown in Fig. 6 and remains very high. 

Before the time t = 0, the disk stays in its axisymmet- 
ric equilibrium state and possesses only an azimuthal motion. 
At t = 0, we switch on the perturbation by adding an asym- 
metric rotating component to the magnetic field. We then let 
the system evolve during more than one thousand orbital rev- 
olutions of the inner edge of the disk. We performed four sets 
of simulations. In the first one, the gravitational potential was 
Newtonian. In the second, we used a pseudo-Newtonian poten- 
tial in order to take the ISCO into account. This is well-suited 
to describing the Schwarzschild space-time. In the third one, 
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Fig. 5. Thin disk approximation in which the ratio H/R never ex- 
ceeds 0.1, on the left for the Newtonian potential and on the right 
for the Schwarzschild geometry. 




Fig. 6. Plasma /3 parameter for the weakly magnetized disk plotted on 
a logarithmic scale. The magnetic field dominates only very close to 
the inner edge, on the left for the Newtonian potential and on the right 
for the Schwarzschild geometry. 

we took the angular momentum of the star into account by in- 
troducing a pseudo-Kerr geometry. And finally in the fourth 
and last set, we performed simulations with a counter rotating 
accretion disk evolving in the Newtonian potential described in 
the first set. 

5.2. Numerical test 

We have tested the pseudo-spectral method used in this paper 
on some standard problems. For instance, we solved the Euler 
equations in cylindrical coordinates for an initial jump in the 
density profile and retrieved the classical results of the shock 
problem with the correct shock speed. However, in order to re- 
duce the oscillations (Gibbs phenomenon) due to discontinuity 
at the shock, we filtered out the high frequency components in 
Fourier-Tchebyshev space. 

In a second set of tests, we launched an m = 1 free wave 
in the magnetized accretion disk, without any perturbation in 
the magnetic field. The wave is slowly damped due to the ab- 
sorbing boundary conditions. The linear analysis elaborated in 
the previous section is in agreement with the non-linear simu- 
lations presented below, at least when the perturbations in the 
variables (p, r, v, B) remain small. For problems without shock 
or discontinuity in the flow, this pseudo-spectral code gives 
very accurate solutions with only a few discretization points. 

5.3. Newtonian potential 

First, we study the behavior of the disk in the Newtonian poten- 
tial. In this case, the Keplerian rotation rate, the vertical and the 
radial epicyclic frequencies for a single particle are all identi- 
cal. However, for a thin gaseous disk, there is a slight difference 
on the order H/R and we have : 

£l k * K r * Kt (94) 

The normalized rotation rate of the star around the z-axis 
is = 0.0043311. Assuming a 1.4M© neutron star, this cor- 



responds to a spin of v, = 100 Hz. The disk's inner boundary 
is located at R\ = 6 and the outer boundary at R2 = 60. The 
orbital angular motion at the inner edge of the disk is Q,„ = 
R x = 0.0680. We normalize the time by dividing it by the 
spin period of the star T„ = jf- = 1450.7. 

The final snapshot of the density perturbation in the disk 
calculated as Ap/po = p/po - 1 is shown in Fig. 7. We clearly 
recognize the spatial separation between the propagating and 
non-propagating regions delimited by the inner Lindblad res- 
onance located at /?" = 23.7. However, the outer Lindblad 
resonance at R°"' = 49.3 is not seen in this run because the 
magnetic perturbation close to the outer boundary is very weak 
compared to its amplitude close to the inner edge. The corota- 
tion resonance, located in the neighborhood of r = 40.0, is also 
not apparent. The reason is that its growth rate is linear and 
very weak; an estimate is found by taking the typical driven 
resonance for a harmonic oscillator as 

d 2 x 

-p- + (1) 2 x- f cos(w t). (95) 

which has a solution given by : 
f 

x{i) = ^—t sin(6>0- (96) 
2 oj 

The numerical value from the simulation is f/2o> ~ 10~ 10 . 
In order to get an appreciable amplitude at the corotation 
point of 10~ 3 of the average density, we would have to wait a 
time At > 1500 T t , which is much more than the total duration 
of our run (about 70 T„). This explains the missing corotation 
resonance at this stage. However, the most interesting features 
are the very pronounced density fluctuations in the innermost 
part of the accretion disk. They correspond to the non-wavelike 
disturbance of the density related to the inhomogeneous solu- 
tion to the Schrodinger type equation (43). They keep their own 
identity even after 1000 revolutions of the innermost orbit and 
rotate at the stellar spin rate. After a few hundred orbital rev- 
olutions, the disk settles down to a new quasi-stationary state 
in which these disturbances persist. This is confirmed by in- 
spection of Fig. 8 in which we have plotted a cross section 
of the density perturbation Sp/po for a given azimuthal angle, 
namely <p — 0. As expected, the density perturbations vanish at 
the disk edges. They behave as predicted by the linear analysis 
in Fig. 3 when allowance is made for the different boundary 
conditions. 

The nonlinearities are weak throughout the simulation. 
Indeed, looking at the Fourier spectrum of the density in Fig. 9, 
the amplitude of each component, plotted vs the mode m on 
a logarithmic scale, is small except for that corresponding to 
the excitation mode. The odd modes are not present. However, 
the weak nonlinearities create a cascade to high even modes 
starting with m = 2. The largest asymmetric expansion coef- 
ficient C m is m = 2, while the next even coefficients roughly 
follow a geometric series with a factor q = 10 -3 , so we can 
write for all m even, C m « q m/2 ~ l C2, until they reach values 
less than 10~ 20 , which can be interpreted as zero from a nu- 
merical point of view. As deviation from the stationary state is 
weak, the amplitudes of these even modes decay compared to 
the previous one, the highest being of course m = 2. Thus, even 
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in the full non-linear simulation, the regime remains quasi- 
linear. Consequently, the parametric resonance phenomenon 
discussed in the previous section is irrelevant at this stage of 
our work. The effects of strong nonlinearity will be studied in 
a forthcoming paper. Note that due to the desaliasing process, 
the modes m > 9 are all set to zero. Note also that the free 
wave solutions leave the computational domain and do not re- 
turn. Only the non-wavelike disturbance produces significant 
changes in the density profile. 




Fig. 7. Final snapshot of the density perturbation Sp/po in the accre- 
tion disk evolving in a quadrupolar perturbed Newtonian potential. 
The disk extends from R l = 6.0 to R 2 = 60.0. The rotation rate of 
the star is Q, = 0.0043311. The time is normalized to the spin pe- 
riod T, = 1450.7. The in = 2 structure emerges in relation to the m = 2 
quadrupolar magnetic perturbation. 
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Fig. 8. Cross section of the density perturbation 6p/po in the 
Newtonian disk at the final time of the simulation. The inner and outer 
Lindblad resonances are shown by vertical lines at r™ = 23.7/49.3. 
Nevertheless the disk feels only the inner resonance, the outer one be- 
ing much too weak because of the very low perturbation amplitude in 
this region. 
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Fig. 9. Amplitude of the Fourier components of the density perturba- 
tion. The odd modes are numerically zero. Due to small nonlinearities, 
the even modes are apparent but with weak amplitude. The compo- 
nents m > 9 are set to zero because of the desaliasing process. 



5.4. Pseudo-Schwarzschild potential 

In order to take into account the modification of the radial 
epicyclic frequency due to the curved space-time around a 
Schwarzschild black hole, we replaced the Newtonian po- 
tential by the logarithmically modified potential proposed by 
Mukhopadhyay (2003). This potential is well-suited to approx- 
imating the angular and epicyclic frequencies in accretion disks 
around a rotating black hole. The radial gravitational field de- 
rived from this potential is given by 



GM, 



1 



R I — R "' s 

1 20 r 



13 r i 

~ " 27 ln Qr-R ms )W, 



,(97) 



Rn 



3 + Z 2 + V(3-Z,)(3+Zi +2Z 2 ) 



Zj = l+(l-a 2 ) 1/3 [(l+a) 1/J + (l-fl) 1/J ] 



1/3 . 



Z 2 = ^3a 2 + z2, 



(98) 
(99) 

(100) 



where a is related to the angular momentum of the star. 

The important feature of this potential is that the radial 
epicyclic frequency for a single particle vanishes at the inner- 
most stable circular orbit (ISCO), which implies that the para- 
metric resonance condition can be fulfilled in the very inner 
part of the accretion disk for different values of the integer n 
but still with a very small radial extension. However, due to the 
small perturbation we remain in the linear regime, and the para- 
metric resonance is only of second order. The disk boundaries, 
the star rotation rate, and the normalization are the same as in 
the previous case. 

The final snapshot of the density perturbation in the disk 
is shown in Fig. 10. The corotation resonance is now located 
at r — 43.7, which differs slightly from the previous simula- 
tion because the orbital velocity is no longer Keplerian due to 
the pseudo-Newtonian force Eq. (97). The non-wavelike distur- 
bance exists between the inner edge of the disk and the inner 
Lindblad resonance. By viewing a cross section, the Lindblad 
radii appear clearly as seen in Fig. 11. Inspection of the Fourier 
coefficients of the density Fig. 12 proves that the linear regime 
is preserved. The level of excitation of the even modes m + 2 
is very weak. 
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Fig. 10. Same as in Fig. 7 but the Newtonian potential is replaced by a 
pseudo-Schwarzschild one. 



Fig. 13. Same as in Fig. 10 but with an azimuthal mode m = 1. 















Fig. 14. Same as in Fig. 11 but with an azimuthal mode m = 1. The 



Fig. 11. Same as in Fig. 8 but for the pseudo-Schwarzschild potential. Lindblad resonances are located at r. 
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12.9/56.0. 



The Lindblad resonances are located at r^' "' = 21.6/45.5. 
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are now excited starting with m = 1 and with decreasing mag- 
nitude. Fig 15. From a qualitative point of view, we can draw 
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Fig. 12. Same as Fig. 9 but for the pseudo-Schwarzschild potential. 
The Fourier coefficients form a decaying geometric series. 
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The results of replacing the m = 2 mode by a more re- 
alistic m = 1 mode, in accordance with the dipolar struc- 
ture of the magnetic field, is shown in Fig. 13. The inner 
Lindblad resonance moves closer to the inner edge and is lo- 
cated at = 12.9. There is less space left for the non-wavelike 
disturbance between the inner edge and the inner Lindblad res- 
onance as can be seen in Fig. 14. Also, all the azimuthal modes 



Fig. 15. Same as Fig. 12 but with an azimuthal mode m 
Fourier coefficients form a decaying geometric series. 



1. The 



the same conclusions as in the case of the Newtonian disk. The 
absence or presence of the ISCO does not affect the physics of 
the resonances, but only has an effect on their amplitude and 
their location. 
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5.5. Pseudo-Ken potential 

We also ran a simulation in which the rotation of the star is 
taken into account. This shifts the location of the ISCO closer 
to the surface of the neutron star. 

We chose a star with an angular momentum of a* = 
0.5. Therefore, the disk inner boundary corresponding to the 
marginally stable circular orbit will be located at /?] = 4.24, 
while the outer boundary is chosen to be at R2 = 42.4. The star 
magnetic perturbation in this run is ten times larger than in the 
previous case, 5b — 0.01. Therefore the computation starts to 
leave the linear regime and reaches a weak non-linear state. 

The corotation resonance at r — 23.8 is not reached, as 
can be seen from Fig. 16. However, the wave and non-wavelike 
disturbances are clearly identified between the inner edge and 
the inner Lindblad resonance. They persist during the whole 
simulation at the low level predicted by the linear analysis. 
The precise location of the Lindblad resonance is determined 
by the cross section view shown in Fig. 17. Looking at the 
Fourier components in Fig. 18, all the even modes are excited 
at a substantial level. The stronger stellar perturbation starts to 
excite modes which are multiples of the injection scale given 
by m — 2. 




Fig. 16. Final snapshot of the density perturbation in the accretion disk 
evolving in a perturbed pseudo-Kerr potential with a = 0.5. The disk 
extends from R[ = 4.24 to R2 = 42.4. The outer Lindblad resonance 
is not on the grid. 



5.6. Retrograde disk 

Finally we checked that the Lindblad resonances disappear for 
a retrograde Newtonian disk. The numerical values are the 
same as in Sec. 5.3, except that the sign of the stellar spin 
is taken to be £2, = -0.0043311 and the magnetic perturba- 
tion 5b = 0.01. 

Looking at Fig. 19, the evolution of the disk does not show 
any Lindblad resonance as expected. The mode m = 2 with two 




Fig. 17. Same as in Fig. 8 but for the pseudo-Kerr potential. The inner 
Lindblad resonance appears clearly at r£ = 19.2 while the outer one 
at rf" = 45.8 lies outside the computational grid 




Fig. 18. Same as Fig. 9 but for the pseudo-Kerr potential. 



spiral arms propagates in the whole computational domain. The 
disk reaches a quasi-stationary state in which the density fluctu- 
ations are well established, Fig. 20. They fit the shape found in 
the linear analysis, Fig. 4. Again, the stronger stellar magnetic 
perturbation affects the high order azimuthal modes, which are 
excited at the same amplitude as the injection mode m = 2, 
Fig. 21. 

Here also we can replace the m = 2 mode by a more real- 
istic m = 1 one, in accordance with the dipolar structure of the 
magnetic field and with, for example, a smaller magnetic per- 
turbation taken again to be 5b = 0.001. The results for the den- 
sity perturbation are shown in Fig. 22. A spiral wave is formed 
and propagates in the whole disk with a rotation rate equal to 
the neutron star spin. Now all the Fourier components are ex- 
cited, the strongest one corresponds to m = 1, and the lighter 
order ones have decreasing amplitudes, as shown in Fig. 24. 

To summarize, in all of the results obtained by these 
runs, the exact structure of the gravitational potential, be it 
Newtonian or pseudo-Newtonian, does not affect the physics 
of the resonances. The main features of these resonances are 
their location very close to the inner boundary of the accretion 
disk. When the applied magnetic asymmetry is small enough, 
the perturbations remain in the linear regime. 

6. CONCLUSION 

In this paper, we studied the effects of an asymmetric rotating 
dipolar or multipolar magnetic field on a magnetized accretion 
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Fig. 19. Final snapshot of the density perturbation in the counterrotat- 
ing accretion disk evolving in a perturbed Newtonian potential. Same 
values as in Fig. 7 apply except for the sign of Q.,. A trailing spiral 
density wave of m = 2 propagates in the entire disk. 
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Fig. 22. Final snapshot of the density perturbation in the counterrotat- 
ing accretion disk evolving in a perturbed Newtonian potential. Same 
values as in the Fig. 7 caption excepted for the sign of fl„. A trailing 
spiral density wave of m = 1 is propagating in the whole disk. 
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Fig. 20. Same as Fig. 8 but for the retrograde disk. The Lindblad res- Fig. 23. Same as Fig. 8 but with a retrograde disk. The Lindblad reso- 



onances are no longer present. 



nances have disappeared. 
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Fig. 21. Same as Fig. 9 but with a retrograde disk. 
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Fig. 24. Same as Fig. 9 but with a retrograde disk. 



disk. This situation is well-suited to LMXBs containing a neu- 
tron star, since they are believed to possess a magnetic field 
that more or less channels the accreting matter to the polar 
caps of the neutron star, depending on its strength. The same 
conclusions can be drawn in the case of a hydrodynamical disk 
evolving in a quadrupolar gravitational potential perturbation. 



Although the flows differ, the consequences of a rotating asym- 
metric field are qualitatively the same. The disk is subject to 
the corotation resonance, the Lindblad driven resonances, and 
if the perturbation is strong enough, the parametric instability. 
We derived an explicit formulation for the regions where the 
blobs should form, given by the resonance conditions Eq. (67)- 
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(69). On the one hand, these blobs can persist during thousands 
of orbital revolutions and can account for the high quality fac- 
tor Q seen in some high frequency-QPOs. On the other, the 
pattern speed of the free wave solutions is closely related to 
the location of the inner boundary of the accretion disk. The 
observed shift in the HFQPOs in response to the accretion rate 
can be explained by these solutions. 

This study should be extended to include viscosity of the 
gas flowing in the disk, which should set a lower limit to the 
magnetic field perturbation amplitude in order to have para- 
metric resonance. In addition, the loss of energy by radiation 
should also damp these oscillations and select only a few of 
them, which will also permit estimation of the light curves of 
such systems. The relation between the simulations presented 
in this paper and observations will be investigated in a forth- 
coming paper, in which we show how to calculate the power 
spectrum density (PSD) of the accretion disk. The prediction 
of the PSD will be crucial for validation of the model. 

We would like to emphasize that this work is only a prelimi- 
nary study undertaken to prove that a rotating asymmetric mag- 
netic field leads to some interesting resonance mechanisms. 
The toy model used here cannot be applied directly to obser- 
vations. Indeed, many aspects, such as the warping/precession 
of the disk expected in 3D simulations, as well as viscosity 
and radiation, should be taken into account in a more realistic 
model. 

Nevertheless, we believe that this idea can also be ap- 
plied to the LMXBs containing a black hole. Indeed, it seems 
that the existence of QPOs in black hole candidates is corre- 
lated with the detection of a jet emanating from these systems, 
(McClintock & Remillard 2003). The Blandford-Znajek pro- 
cess (Blandford & Znajek 1977) converting the rotational en- 
ergy of the black hole into the launching of the jet should also 
be responsible for a connection between the accretion disk and 
the black hole via the magnetic field. Therefore, the model pre- 
sented in this paper has the potential to explain the QPOs seen 
in the BHCs (Wang et al. 2003). 

Acknowledgements. This research was carried out in an FOM projec- 
truimte on 'Magnetoseismology of accretion disks', a collaborative 
project between R. Keppens (FOM Institute Rijnhuizen, Nieuwegein) 
and N. Langer (Astronomical Institute Utrecht). This work is part 
of the research programme of the 'Stichting voor Fundamenteel 
Onderzoek der Materie (FOM)', which is financially supported by the 
'Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)' . 

I am grateful to the referee Lev Titarchuk for his valuable com- 
ments and remarks as well as to Jean Heyvaerts and John Kirk for 
carefully reading the manuscript. 

References 

Alpar, M. A., & Shaham, J., 1985, Nature, 316, 239 

Anzer, U., & Borner, G., 1980, A&A, 83, 133 

Anzer, U., & Borner, G., 1983, A&A, 122, 73 

Blandford, R. D., & Znajek, R. L., 1977, MNRAS, 179, 433 

Belloni, T, Psaltis, D., & van der Klis, M., 2002, ApJ, 572, 392 

Canuto, C, Hussaini, M., Quarteroni, A., & Zang, T, Spectral 

Methods in Fluid Dynamics , Springer Verlag (Berlin and New 

York), 1990. 



Frieman, E., & Rotenberg, M., 1960, Reviews of Modern Physics, 32, 
898 

Ghosh, P., Pethick, C. J., & Lamb, F. K., 1977, ApJ, 217, 578 

Ghosh, P., & Lamb, F. K., 1978, ApJL, 223, L83 

Ghosh, P., & Lamb, F. K., 1979a, ApJ, 232, 259 

Ghosh, P., & Lamb, F. K., 1979b, ApJ, 234, 296 

Goldreich, P., Tremaine, S., 1979, ApJ, 233, 857 

Gottlieb, D., & Orszag, S., 1977, Numerical Analysis of Spectral 

Methods: Theory and Applications SI AM 
Kato, S., 2001a, PAS J, 53, 1 
Kato, S., 2001b, PASJ, 53, L37 
Kato, S., 2004, PASJ, 56, 905 

Kluzniak, W., Abramowicz, M. A., Kato, S., et al., 2004, ApJL, 603, 
L89 

Lai, D., 1999, ApJ, 524, 1030 

Landau, L., & Lifshitz, E., 1982, Cours de physique theorique, Tome I 

mecanique, Editions Mir Moscou. 
Lee, W. H., Abramowicz, M. A., & Kluzniak, W„ 2004, ApJ, 603, L93 
Lindblad, P. O., 1974, IAU Symp. 58, 399 
Mauche, C. W., 2002, ApJ, 580, 423 

McClintock, J. E., & Remillard, R. A., 2003, astroph/0306213 
Mauche, C. W., 2002, ApJ, 580, 423 

Miller, M. C, & Lamb, F. K, & Psaltis, D„ 1998, ApJ, 508, 791 
Morse, & Feshbach, 1953, Methods of Theoretical Physics, New- 
York: McGraw Hill 
Mukhopadhyay, B., & Misra, R., 2003, ApJ, 582, 347 
Psaltis, D., Belloni, T, & van der Klis, M., 1999, ApJ, 520, 262 
Schnittman, J. D., & Bertschinger, E., 2004, ApJ, 606, 1098 
Shirakawa, A., & Lai, D., 2002a, ApJ, 564, 361 
Shirakawa, A., & Lai, D., 2002b, ApJ, 565, 1 134 
Smirnov, V., 1989, Cours de mathematiques superieures, Tome 3 

vol. 2, Editions Mir Moscou. 
Stella, L., & Vietri, M., 1998, ApIL, 492, L59 
Stella, L., & Vietri, M., 1999, Physical Review Letters, 82, 17 
Stella, L., Vietri, M., & Morsink, S., 1999, ApJ, 524, L63 
Titarchuk, L., & Fiorito, R., 2004, ApJ, 612, 988 
Titarchuk, L., & Wood, K, 2002, ApIL, 577, L23 
Titarchuk, L., Lapidus, I., & Muslimov, A., 1998, ApJ, 499, 315 
Titarchuk, L., 2003, ApJ, 591, 354 
van der Klis, M., 2000, ARA&A, 38, 717 
Vietri, M., & Stella, L., 1998, ApJ, 503, 350 
Wang, D., Ma, R., Lei, W. & Yao, G„ 2003, MNRAS, 344, 473 
Wijnands, R., 2001, Advances in Space Research, 28, 469 
Zhang, W., Smale, A. P., Strohmayer, T. E., & Swank, I. H., 1998, 
ApJ, 500, L 171 

Appendix A: Eigenvalue problem for the 
Lagrangian displacement 

A.1. The Lagrangian displacement £ 

The linearization of the MHD equations is done by introduc- 
ing the Lagrangian displacement vector which is related to 
the Eulerian and Lagrangian perturbation of the fluid physical 
quantities, such as the velocity, the magnetic field, and the den- 
sity. These perturbations are denoted respectively by 6 and A. 
By definition, the Lagrangian perturbation of the velocity is 

Av=v-v = ^ = ^+v-V£ (A.l) 

Moreover, the Eulerian perturbation is related to the 
Lagrangian one by 

tf = A-£-V (A.2) 
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Therefore, the Eulerian perturbation of the velocity reads 



6v=vi = ^+vVf-f Vv 
ot 



(A.3) 



The boundary conditions at the disk- vacuum interface are 



B 2 

<p + —> = 
2po 
n B = 

iiA < B > = 



(A.4) 

(A.5) 
(A.6) 



where < X >= X vac - Xdisk denotes the jump at the interface 
and n is a unit vector normal to the plasma interface. 

Let's start by perturbing the mass conservation equation by 
using the equilibrium condition div (poVo) = 



-p- +v ■ Vpi +pidivv = -div(p Vi) 
Ot 

-divU^- +Po(v- V£-£- V 



(A.7) 



From the tensorial calculus formalism, the right hand side can 
be transformed after tedious algebra into 



+ v • Vpi + pi div v = |- (-div (p £)) 
+v • V (-div (p £)) - div (p f) div v 



(A.8) 



Because the Lagrangian displacement is arbitrary, we make the 
following identification 

Pl =-div(po£) (A.9) 
To continue by perturbing the induction equation, we get 
ggi 



+ v • VB] + Bi divv - Bi • Vv = 
B • Vvi - Bq divvi - V\ ■ VB 



We transform each term on the right hand side from the equilib- 
rium state rot (v A B ) = 0. The unperturbed magnetic field is 
divergenceless neither is the perturbed one, div Bo = divBi = 
0. Thus, we obtain 



Bo • Vvi = -(Bo • V0 + v • V(B„ • V£)+ 
ot 

div v B • V£ - div (B £ • Vv ) 

-Bodivvi = -— (B div£) - v • V(B div£)+ 
Ot 

(v • VB )div£ + B £Vdiv v 
d 

-v, ■ VB = --(£ • VB ) - v ■ V£ ■ VB + 

f • V [vo • VB ] 

By introducing the vector 
Q = rot A Bq), 



(All) 



(A. 12) 



(A. 13) 



(A. 14) 



the sum of the three terms can be cast after tedious tensorial 
calculus into the form 



^ + v • Vg + gdivvo - Q ■ Vv 
ot 



So we make the identification B\ - Q. This is the perturbation 
of the magnetic field induced by motion in the disk and does 
not include the magnetic field perturbation emanating from the 
star. 

Finally for the pressure perturbation, we develop the adia- 
batic condition to first order. We have 

^P- +vo • Vpi +ypidivv = -yp divvi -V\ • Vp (A.16) 
ot 

Introducing n = -y po div £ - £ • V po, the right hand side is 
dn 



dt 



+ Vo • V7T + y n div Vo 



and thus the pressure perturbation of the gas is 
Pl =5p= -yp div£ Vp 



(A. 17) 



(A. 18) 



Equations (A.3), (A. 14), (A. 18) are the perturbed quantities in- 
duced by the motion in the fluid. However, the magnetic field 
perturbation possesses another variation that is not imposed by 
the fluid but is introduced independently by the neutron star. 
This known perturbation 5B» has to be added to the momentum 
equation. In the realistic case of a dipolar magnetic field, we 
should expect rot B« = 0, as well as rot SB* = 0. Nevertheless, 
in our 2D numerical simulations, we do not keep this property 
for the stellar magnetic field because in this particular geome- 
try, we impose straight field lines aligned along the rotational z- 
axis. Therefore, in order to treat the problem in a self-consistent 
manner we keep this non-vanishing curl term in the linear anal- 
ysis already in order to facilitate the connection linear analysis- 
numerical simulations. To get the eigenvalues system satisfied 
by £, we put expressions (A. 14, A. 1 8) into the momentum equa- 
tion (16). 



(A. 10) A.2.PDEfor£ 



Replacing all the perturbed quantities by their expression in- 
volving the Lagrangian displacement and keeping only the first 
order term, satisfies a second order linear PDE given by 

dv\ 

Po ^- + Po vo • Vvi + po vi • Vv + pi v • Vv = 
ot 

— (rotB Q A (Bi + SB*) + rotB x A (B + 

po 

SB,) + rot SB t A (B + B x + SB,)) (A. 19) 

On the left hand side, we have for each term 

^ =Po0 + Pov o -vf -p f -Vvo (A.20) 

po vo • Vvi = po v • V— + po v • V(v • Vf - f • Vv ) 
ot 



(A. 15) 



Po vi • Vv = po — • Vv + po (v • V£ • V) v 
Ot 

-po (f • Vv • V) v 
Pi v • Vv = -div (p £) v • Vv 



(A.21) 

(A.22) 
(A.23) 
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Finally, summing these four terms gives 

<9 V 1 TT „ _ 

Po-^r + Po vo • vvi + po Vi • vv + pi v • Vv = 
at 

Po -gp + 2p v • V— + div [p v (v • V)£ - p f (v ■ V)v ] 



(A.24) 



Recalling that g — B\ and p\ — n, we rewrite 



roiBo AQ + rotQAB () = -V(B g) + Bo Vg + g VB 

(A.25) 

Introducing the opposite of the total pressure variation, 
gaseous and magnetic, by II = y po div£+£- Vpo - ^ BqQ, the 
Lagrangian displacement for the fluid motion in the accretion 
disk satisfies the following second order linear partial differen- 
tial equation 

Po^ + 2pov -V- + 
div [p v (v • V)f - po £(vo ■ V)v ] - 

vn - — Bo • Vg - — g • VB + 

Ho Ho 

div (po £)(go + 6g) -poSg = 

— [rot (B + Q + SB,) A 6B,+ 
Po 

rotSB„ A(B +g)] (A.26) 

In the expression above, we allow for a perturbation in the grav- 
itational field induced by the star and denoted by Sg. However, 
in the following analysis, we remove this freedom and only 
keep a perturbation in the magnetic field. 

A. 3. Projection onto the cylindrical coordinate system 

In the stationary state, the disk orbits around the compact object 
without mass in- or outflow from the edges and without preces- 
sion. From now on, we drop the subscript for the stationary 
quantities, which should not lead to any confusion. Thus the 
velocity in the axisymmetric stationary state is simply 



v = v v e v = rClep 



We find the following expressions 
(v • V)v 



-r Q 2 c r 



(v-V)£ = Q 



dtp 



(A.27) 



(A.28) 



(A.29) 
(A.30) 



(£-V> = -n^e r +£-V(rn)e v 
We also give the explicit formulas for each term in Eq. (A.26). 

a 2 £ d% d% 



dt 1 



-^e z (A.31) 



d£ d 
2pvV-|=2pQ- 



dip 



dip e * dtp* 1 



(A.32) 



„„ an l dn du 

Vn = — e r + - — e v + — e z 

dr r dip dz 

-div(p£(v- V)v) = 
[rQ 2 div(p^)+p^V(rn 2 )] c r +p^0 2 ^ 

div(pv(v-V)£)=pQ 2 ff 



dip 2 



= pfl 2 



iw- & - 2 - ]e -' 



dp 2 ^ + 2 &p~j e * + ~dy 2 ~ ez 

B Vg = 

9Qr] 



dip 

d 2 % z 



B r 



dQ r Bp dQ r B v Qp 



dr 


r dip 


dQ v + 
dr 


B v dQp 


r dip 


+ 


Br ^r~- 





+ B 



dz 



+ B 



z dz 

r dip z dz 

Q VB = 

dB r Qp dB r QpBp ^ dB r 
dr~ 



r dz 



r dip r 
n dBp Qp dBp QpB r 

Qr + — TT- + — " 

dr r dip r 



+ Qz 



dBp 
~~dz 

dB 7 



— . '^L — . Q 

dr r dip z dz 



pSg=p (5g r e r + Sgp ep + 5g z e z ) 
-div (p£)(g + Sg) = 
-div (p f) [(g r + Sg r ) e r + Sgp ep + (g z + Sg z ) e z ] 
rot (B + Q + SB*) A SB t = 



SBl \-(B r + Q r + SBl) ■ 
j- (B z + Q z + SBl) 



SB? { d i w x 

-f \j r {r(Bp + Qp + SB?))- 



dp 



(B r + Q r + SBl) 



6 -^{l{r(Bp + Qv + SB?)). 



— (B r + Q r + SBl) 

dtp 
1 d 

SB : A-—(B Z + Q Z + SBD- 



d_ 

dz 



(Bp + Qp+ SB?) 



1 d 

SB?\-—(B Z + Q Z +SBD- 
[rdtp 

j z (Bp + Qp+SB?)^- 
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(A.33) 

(A.34) 
(A.35) 



(A.36) 



(A.37) 
(A.38) 

(A.39) 
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5B[ \ -{B r + Q r + 5B r ,)- 



— (B z + Q Z + 6BI) 
or 



(A.40) 



B, r 



(B z + Qd 
Q v 



rot 5B« A(B + Q) 

d6B r , 36 Bl' 



I85B: d6Bl\ 
\ dz dr I 
36B 



r 

B r + Q, 



Ur6Bt) 
or 



dtp 

Or dtp 



.'I 88BI 
(B z + Q Z )\- 



85 BV 



(B r + Qr) 



dtp 


dz /. 


1 d5B% 


dSBl 


r dtp 


dz 


iosb: 


98Bl\ 


\ dz 


dr ) 



(A.41) 



The perturbations in the disk magnetic field induced by the 
fluid are in cylindrical coordinates 



Qr 



Id, x d 

■ B<p - B r J - — (ft B r ~ ft B z ) 



r dtp 



Qip = ~q z (ft B z - ft By) - — (ft By - ft B r ) 

1 d 1 d 

Q - z- ( r (6 B r - ft B z )) - z —(g v B z - ft B v ) 



(A.42) 
(A.43) 



dr 



r dtp 



Projecting on the 3 axis, and taking the stationary condition 
into account 



B 1 



(A.45) 



The radial component of the Lagrangian displacement satisfies 

D 2 & „ „ £>ft dU 1 dp ,. , „ 
Dt z Dt dr p dr 



prft V(Q 2 )-- 



B ^Q L + B 1 dQ L 
r dr r dtp 



BpQp dQ r dB r Q v dB r dB r 

Z + £>z + Qr -7— H 7T- + Hz "T - 



dz 



dr 



dtp 



dz 



(p-div(p€))5g r + — 



1 

j"0 



6Bl\-(B r + Q r + 6B0- 



-(B Z + Q Z + 6B: 
dr 



8Bt (d , „ x 



— (B r +Q r + 6B 
dtp 



+ 



(b z + &) 



«9z 



<9r 



The azimuthal component is given by 



Dft i dn 



l 



O r — 1 1 h 

or r dip r 

z <9z 



Q v dB v B r 

Qr I ^ 



dr 
1 

y"0 



(p-div(p£;))Sg lf ,+ 
^{|(r(B, + G,+«0)- 

^-(B r +Q r + 5s;)J - 
Br + Qr (UrtBty 



dr 



(B z + Qz)\~ 



r dtp 



d5B[ 
dtp 
dSBt 

~ ~d~r 



(AM) 



and finally the vertical component satisfies 



(A.44) ^ rPj^ _dn _j_ 



dQ z B v dQ z dQ z 

Dr — 1 " r D, — h 

dr r dtp dz 
„ dB z dB z „ dB z 

Qr + + Q? 

dr r dtp dz 
(p-div(p€))6g z -div(pi;)g z + 



1 

y"0 



ld_ 

r dtp 



(B z + Q Z + 6BD- 



-(B^ + Q^+SBt) 



SB* 



SBl\ — (Br + Qr + SS0- 



-(B z + Q z + SBl)\ + 



. 1 dSBl dSBt 
(B v + Q v ) - 



(Br + Qr) 



r dtp 
d6B 



dz 

[ dz dr J 



(A.48) 



Eq. (A.46)-(A.47)-(A.48) are the fundamental equations 
governing the perturbed motion in the disk due to a small grav- 
itational and magnetic field. 

A. 4. Orbital and vertical decoupled oscillations 

To elucidate the meaning of these equations, we focus on the 
(A.46) disk response when motions in vertical direction and orbital 
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plane are decoupled. Let's start with the Lagrangian displace- 
ment contained in the orbital plane of the disk by setting J; = 
and £. = 0. We also assume that perturbations are axisymmet- 
ric in order to avoid meaningless complications so that ^ = 0, 
corresponding to m = modes. Then the magnetic perturba- 
tions read 



Qr = 



(A.49) 
(A.50) 

(A.51) 



and the opposite of the total pressure perturbations 

B z Q z 



1 d dp 

n = y/»-^-(rf r )+?r/ 
r or or 



Po 



Keeping only the leading terms in Eq. (A.47) we find a simple 
relation between radial and azimuthal displacement, namely 



^ + 2Q^ 
Dt 2 Dt 







(A.53) 



This can be integrated directly. Dropping the constant of inte- 
gration we find 



Dt Sr 



(A.54) 



After straightforward algebra and making use of the equi- 
librium condition Eq.(6) and the relation (A.54), the radial 
Lagrangian displacement (A. 46) satisfies the following second 
order partial differential equation 



D% 
Dt 2 



d_ 

dr 



B z 6Bl\ 1 3 



Po 



, 3 (SB? />• 8B, 

/J.Q 



dr \ 2p 



d 18B% 8B Z 
dr\ pio dr 



(A.55) 



k 2 = 4 Q 2 + r is the radial epicyclic frequency of the os- 
cillations. The sound speed and the Alfven velocity associated 
with the vertical magnetic field are respectively 



7P_ 
P 

Pop 



(A.56) 
(A.57) 



In the vertical direction, we neglect the variation in the (r, if) 
plane implying that ^ = ^ = 0. We set also (^ r ,^) = such 
that the magnetic perturbations are 



Q r = -— (B r £) 
oz 



and the opposite of the total pressure perturbations 

df z dp B r Q r 

n = rp-7-+fz-7 

oz oz no 



(A.58) 

(A.59) 
(A.60) 



(A.61) 



After straightforward algebra and making use of the equi- 
librium condition Eq.(7), the vertical Lagrangian displace- 
ment (A. 48) satisfies the following second order partial differ- 
ential equation 



D 2 & 
Dt 2 

d_ 

d~z 



d_ 

dz 



, 2 2 N B r 6B r A d& 
p{c]+c 1 ar ) + '-]-£■ 



Po 



dz 



+ P«Uz = 



'SB^+SBf B r 6B: 



2^o 



po 



+ d_l5W 1 dB L \ 

dz \ Ho dz 7 



(A.62) 



The Alfven velocity associated with the radial magnetic 
field is 



(A.52) -2 = B r 

ar POP 



(A.63) 



A. 5. Derivation of the inner boundary condition 
Eq. (51) 

In the numerical runs presented in this paper, we use non re- 
flecting boundary conditions which leads to a vanishing density 
perturbation at the inner edge of the disk, namely dp - 0. In or- 
der to compare the linear analysis with simulations, we adopt 
the same conditions in the both cases. Neglecting the displace- 
ment in the vertical direction, density perturbations are related 
to the Lagrangian displacement by 



5p = -div(p£) = -- ^-(rptr)-i-pt<p 
r dr r 



(A.64) 



To the lowest approximation given by Eq. (A.54), the radial and 
azimuthal displacement expanded by Eq. (32) are related by 



&*-2i-£ 



(A.65) 



Injecting into the density perturbation, we have at the inner 
edge of the disk r - R\ 



d Q. - 

Sp= —(rp£r) + 2m—p£r = 
dr cj 

We recast the terms in such a way to obtain 

dh ( d Q \ - 

dr \dr ro> 



(A.66) 



(A.67) 



evaluated at r — R\. This is the inner boundary condition given 
byEq.(51). 

Appendix B: Approximate solution to the 
Schrodinger type equation 

In this appendix, we recall a method to find approximate solu- 
tions to the following Schrodinger type equation with a source 
term f(x) 



y"(x)+p(x)y(x) = f(x) 



(B.l) 



The functions / and p are continuous in the interval [a,b], 
which includes the origin x = 0, (a < 0, b > 0). Moreover, 
p has one pole in this interval corresponding to the point x — 
where it changes sign, p(0) = and p'(0) + 0. 



(B.13) 
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B.1. Homogeneous solution B.1 .2. p(x) < 

We start to look for solutions to the homogeneous part of Then<jJi(x) > 0. The solution of Eq. (B.ll)is 
Eq. (B.l) when f(x) = 0. The reference equation is Airy's ^ 2 /3 

equation given by wi (x) = ± - J yj-p(s)ds 

w (x) xw(x) - (B.2) minus sign applies if p(x) < for x < 0, otherwise the 

We then look for approximate solutions to Eq.(B.l) which plus sign applies. 

can be cast into the following form, see for example Note that the function o>, is strictly monotonic. co[ never 
Smirnov (1989) changes sign and never vanishes. Consequently, the ampli- 

tude A given by Eq. (B.10) remains finite in the entire inter- 
y(x) = A(x)w(coi(x)) (B.3) V al[a,fo]. 

w is an exact solution of Airy's equation. The functions A 

and (±>\ are chosen such as to satisfy the homogeneous part B.1 .3. Conclusion 
of Eq.(B.l) within a prescribed accuracy. The first and second 

, . .. t , , Finally, two linearly independent solutions of the homogeneous 

derivatives of y(x) are •" 1 r b 

part of Eq. (B.l) can be expressed with help on the Airy func- 

y'(x) = A'(x)w(oj 1 (x)) + A(x)a>[(x)w'(a)i(x)) (B.4) tions of the first and the second kind, respectively Ai and Bi 

y"(x) = A"(x)w(coi(x)) + 2A'(x)a)[(x)w'(o)i(x))+ 1 

1 yi(x) = — = Ai(o)i(x)) (B.14) 

+ A(x) [co'{(x)w'(coi(x))(a>' 1 ) 2 (x)w"(coi(x))) (B.5) VKWI 

However, due to Eq. (B.2), w"(a>i(x)) = a>i(x) w{u>\{x)) and yi{x) = ■ Bi{o)\(x)) (B.15) 
therefore the second order derivative can be transform into V^i^l 

„ , .// . , /s? , . s ,~ ., , //s , ss.„ The associated Wronskian is given by 
y" = (A"+A(w' 1 ) 2 wi)w(wi) + (2A'w' 1 +Aw' 1 ')w'(wi(x))(B.6) 5 3 

t t sign(w'j) 

Inserting this expression into Eq. (B.l), (remember that we W(yuy 2 ) =yiy 2 -y'iy2 = (B.16) 

look for homogeneous solutions with f(x) = 0), we get T ., . , . . ,. 

It never vanish, proving that the two solutions remain linearly 

(A"+A (<jj\) 2 (jJi+pA) w(a»i)+(2 A' u\+A co") w'(o>i(x)) = 0(B.7) independent in the entire interval [a, b]. 
The Airy functions w and W being linearly independent, A and 

coi satisfy the following system B - 2 - Inhomogeneous solutions 

A" + A (to' ) 2 lo\ + p A = (B 8) Knowing the general solutions to the homogeneous part, a pe- 

, , „ culiar solution to the full Eq. (B.l) with source term is given by 

2 A <u, +Aoj { =0 (B.9) theformula 

Assuming that the amplitude A is positive, Eq. (B. 9) is solved C x y\ (x)y 2 (t) - ydOyiix) 

by ' ' Mx) -J wi^W) 



f(t)dt (B.l 7) 



K Consequently, the most general approximate analytical solu- 

~ ^\(d' x (x)\ (B. 10) t - on tQ Qur m j t j a j p ro blem Eq. (B. 1) to the lowest order is 

where K is constant. Neglecting the second derivative A" x y( x ) = c \y\{x) + C 2 y 2 {x)+ 

which is one order of magnitude less than the other terms in s ig n ( w ' 1 ) 7 r C ( yi (x)y 2 (t) - yi(t)y 2 (x))f(t)dt (B.18) 
Eq. (B. 8), it reduces to J 



Eq. (B.8), it reduces to 

(0j\) 2 Uy + P = (B.ll) 



As usual, the constants C\,C 2 are defined by properly fitting the 
imposed boundary conditions at x = a and x = b. Applied to 
This relation implies that the function toi has to satisfy the the linearized MHD accretion disk, we look for solutions which 
property p(x) a)i(x) < 0. We have to distinguish the following remains bounded for large radii. This gives the first bound- 
two cases. ar Y condition. The second boundary condition is derived by an 

appropriate density or pressure perturbation at the inner edge. 
(See end of the appendix) 

B.1.1.p(jc)>0 

It follows that u)i(x) < 0. The solution of Eq. (B.ll) is Appendix C: Derivation of the resonance 

conditions 



wi(x) 



3 r x 

: 2 Jo 



p(s) ds 

o 



2/3 

(B.12) 



C.1. Corotation resonance 



The plus sign applies if p(x) > for x > 0, otherwise the minus The corotation resonance is defined by the radius r c 
sign applies. where co t (r c ) = 0. Actually, this equation possesses two so- 
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lutions corresponding to oj(r c ) = + mc ^\ The width of this re- 
gion is of the order of the disk height 0(H). For very thin disks, 
this discrepancy can be neglected and the two solutions merge 
together in an unique corotation radius given by o>(r c ) = 0. In 
other words, we assume in this case that oj = to,. However, in 
our numerical application, the separation between the two coro- 
tation radii is large enough to be resolved. Nevertheless, in the 
derivation of the corotation resonance below, we will assume 
that corotation is achieved when a>(r c ) = m (Q, - £2(r c )) = 0. 
For the detailed study of both corotation, it can be shown that 
the leading term in the equation satisfied by the radial displace- 
ment Eq. (25) is given by 



r 



1 + 



d ^(rpc 2 maz ) m 2 c 2 



2 2 
m Cmaz 



din r 



+ ■ 



or 



dr 2 ' 
dr ' 



-AD. 



oj- \ m ci, n , d o 

r or 



m 2 c 2 maz B z 5Bl d 



pr 1 - 



(CI) 



Indeed, it is obtained from Eq. (25) by replacing the convective 
derivative D/Dt by -ioj, % v by Eq. (37) and neglecting sec- 
ond order terms. Keeping only the leading divergent terms in 
the coefficients of the ordinary differential equation (C.l), we 
obtain 

m 2 c maz g, tf) + ™ 2 c maz £_ I J_ 

r 2 col r r 2 dr\ oj 2 



z mc maz^ d I 1 | ,- ( „ ) _ 

r dr\oj 2 ." ] 



™ 2 <Zaz 

pr 2 



B 7 5Bl d I 1 



dr \ojI 



We introduce the new independent variable 
r - r r 



x = 



(C2) 



(C3) 



Expanding w, to the first order around the corotation radius r c 
we have 



Oi,{f) = 0)*(r c ) + (r - r c ) oi'„{r c ) + o{r - r c ) = 
x r c Oj' t (r c ) + o(x) « a x 

To this approximation, we have to solve 

^(i)--^(x)-4-^if r (i) = 



(C4) 



-2- 



P C maz " 

This is of the form 



B z (r c )6Bl(r c ) 

y"0 



y"{x)-- x y{x)- b -y{x)=- x 



(C5) 



(C6) 



with£ = 4^ 

111 X 

and d = -2—1 



> because ojx- m(Q, t - fl) (r- r c )/r c > 
- — B ' < - rc '> SB < ( ^\ Making the change of vari- 



able t = 2 ybx and introducing the new unknown v(f) by y(t) 
/ 3 v(t), v(t) satisfies the modified Bessel equation of order 3 



v»(,) + I v '(r)_(l + i)v(0 = 



This is solved by 

v(t) = C 1 I 3 (t) + C 2 K 3 (t) 



(C7) 



(C8) 



Thus, the complete most general solution to Eq. (C.5) for which 
a particular solution is easily found to be a constant equal 



£.(*) = C] x 3 ' 2 7 3 (2 yfbx) + C 2 x V2 K 3 (2 Vfojc) 

m B z (r c )5Bl(r c ) 



2Qwp 



(C9) 



Finally, near the corotation radius, the Lagrangian displace- 
ment which remains bounded needs C\ = 



£.(*) = C 2 x V2 K 3 (2 yfbx) 



m B z (r c )SBl(r c ) 



2Q,ojp 



(CIO) 



The density disturbance induced in the disk by the rotating 
gravitational perturbation is then to the lowest order 



Sp 
P 



div(p£) Id m 

= -i-(rp&) j 

p pr dr ro)% 



m B z (r c )6Bl(r c ) 



rp 



y"0 



(C.ll) 



The displacement Eq. (CIO) is continuous and differen- 
tiable everywhere. Thus, the first term on the right hand 
side has a finite value. The second term on the RHS 
needs a special treatment. Indeed, when r approaches r c the 
numerator(-— B ' ( - rc>6B ' (rc) + 2Q.o>£ r ) and the denominator oj 2 
vanish as well, leading to an undetermined expression of the 
form 0/0. To find the behavior near r c we introduce the func- 
tion f(r) = —— Bz(r) 6B ' (r \ We note that near the corotation, 
(f(r) + 2D co^r) behaves as f(r) - f(r c ) = f'(r c )(r - r c ) 
with f'(r c ) + 0. Thus we conclude that it approaches zero as x. 
Therefore the density perturbations in the disk diverge as 



dp 
P 



mf'(r c ) 
oj' t (r c ) 2 r 2 x 



(C12) 



The divergent term in the density perturbation Eq. (C. 12) tends 
to infinity as i. This result is consistent with the conclusions 
drawn by Goldreich & Tremaine (1979) for a hydrodynamic 
disk without self-gravity. We just need to replace the sound 
speed c 2 by the fast magnetoacoustic wave speed c^ flz and the 
gaseous pressure by the total (gaseous+magnetic) pressure. 

C.2. Lindblad and parametric resonance 

By choosing an appropriate origin of time to, the generalized 
Mathieu equation (66) can be recast in the more convenient 
form as follows 



f {t) + [ K 2 +h cos( r t)\ f(f) = / cos(yf) 



(C.13) 
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where we have introduced the comoving excitation fre- 
quency y = Q. - Q. t . It is a second order linear ordinary dif- 
ferential inhomogeneous equation for the variable The ho- 
mogeneous part corresponds exactly to Mathieu equation. The 
related solutions known as Mathieu functions are well stud- 
ied in mathematical physics textbooks, see for example Morse 
& Feshbach (1953). These authors show that the solution be- 
comes unbounded when the excitation frequency y lies in a 
close neighborhood to a multiple of twice the harmonic fre- 
quency k, namely n\y\ ~ 2 k where n is an integer. The size 
of the neighborhood depends on the amplitude of the modu- 
lation h. Because the largest growth rate is reached when the 
relations holds exactly, we only keep this fastest exponential 
growing solution. This gives the parametric resonance condi- 
tion as exposed in Eq. (69). 

When the amplitude of the modulation is weak compared 
to the driven force h <sc /, the parametric growth rate can be 
neglected. Indeed, to this approximation, Eq.(C.13) reduce to 
the well known driven harmonic oscillator 

r(0 + « 2 f(0=/cos(y0 (C.14) 

Resonance occurs when \y\ = k. In that case the most general 
solution reads 

ft 

£(r) = A cos(Kt) + B sm{Kt) + — sin(Kf) (C.15) 

While in the parametric resonance the growth is exponential, 
here it is only linear with time. Therefore, in the early stage of 
the full solution of Eq. (C.13), this latter form is a good approx- 
imation. On a very long timescale, it converges asymptotically 
the Mathieu function with exponential growth rate. 

Note that the driven resonance (called Lindblad resonance 
in the main text to refer to the same physical process arising in 
galactic dynamic and described by Lindblad (1974)) is a spe- 
cial case of the parametric resonance for n = 2. As a conse- 
quence, we could enclose the driven resonance condition into 
the parametric one. Nevertheless, in order to keep track of the 
difference between linear and exponential growth rate, this dis- 
tinction will remain throughout the paper. Moreover, in the nu- 
merical application shown in the simulations, the parametric 
resonance does not appear because of the weak magnetic field 
disturbance. The timescale of the growing parametric mode ex- 
cess by several orders of magnitude the duration of the simula- 
tions. Only the Lindblad resonance is relevant in this case. 



